GridLib.f90 Source File
manage input and output of grids
Source Code
!! manage input and output of grids !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL ! !### History ! ! current version 1.9 -18th March 2025 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 0.1 | 02/Nov/2007 | Original code | ! | 0.2 | 30/Nov/2008 | extract grid of specified date | ! | 0.3 | 10/Dec/2008 | export grid to netcdf dataset | ! | 0.4 | 13/Dec/2008 | Add support to read ESRI GRID | ! | 0.5 | 18/May/2009 | Add support to grid_mapping (makes use of GeoLib) | ! | 0.6 | 03/Jul/2009 | export to file (ESRI GRID) | ! | 0.7 | 07/Apr/2010 | added GetDtGrid | ! | 0.8 | 19/May/2011 | increased decimal digits in header of exported esri grid | ! | 0.9 | 17/Jun/2011 | add support to read netcdf file from WRF | ! | 1.0 | 11/Mar/2012 | added GetFirstDate and GetTimeSteps | ! | 1.1 | 17/Jun/2013 | added support for grid with non regularly spaced coordinates in netcdf dataset | ! | 1.2 | 15/Dec/2016 | GetXYSizesFromFile modified to detect x and y coordinates by units | ! | 1.3 | 18/Feb/2019 | GetXYSizesFromFile modified to detect coordinate convention lat-lon or lon-lat | ! | 1.4 | 20/Apr/2021 | ScanTimeStringAsString to parse reference time in netcdf dataset | ! | 1.5 | 11/Apr/2022 | SyncTime to return lower bounding time step of a given datetime ! | 1.6 | 06/Apr/2023 | comments reformatted to adhere FORD documenter standard | ! | 1.7 | 08/Nov/2024 | increased number of decimal places of header in exported ESRI ASCII grid | ! | 1.8 | 21/Jan/2025 | Easting and Northing used to find x and y coordinate | ! | 1.9 | 18/Mar/2025 | minor bug fixing in reading coordinate reference system in netcdf files | ! ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! ! This file is part of ! ! MOSAICO -- MOdular library for raSter bAsed hydrologIcal appliCatiOn. ! ! Copyright (C) 2011 Giovanni Ravazzani ! !### Module Description: ! set of fortran routines to manage input and output of grids. ! The module supports different file formats: ! _NetCDF_ raster layer CF 1.0 compliant (Climate and Forecast metadata standard) ! _ESRI ASCII_ grid ! _ESRI BINARY_ grid ! ! The module introduces two new variable types: ! `grid_real` to store floating point data, and ! `grid_integer` to store integer data. ! ! Spatial conventions in netCDF file: !``` ! y ! yDim ^ N ! | ! | W E ! | ! | S ! 1 |---------> x ! 1 xDim !``` ! !index conventions: GRID(x,y) ! ! Spatial conventions used in GridLib: ! !``` ! 1 jdim ! 1 |---------> j ! | N ! | ! | W E ! | ! | S ! idim v ! i !``` ! !index conventions: GRID(i,j) ! To initialize a new grid variable, either integer or float, you can: ! ! read grid from NetCdf file passing the standard name: !```fortran ! TYPE(grid_real) :: new_grid ! CALL NewGrid (new_grid, NetCdfFileName, stdName='air_temperature') !``` ! ! read grid from NetCdf file passing name of variable: !```fortran ! TYPE(grid_real) :: new_grid ! CALL NewGrid (new_grid, NetCdfFileName, variable = 'T2m') !``` ! ! If NetCdf archive is multitemporal, that is contains more grids, ! one for each time, you can retrieve the grid for a given time: !```fortran ! TYPE(grid_real) :: new_grid ! TYPE (DateTime) :: ncTime ! ncTime = '2007-03-05T11:00:00+00:00' ! CALL NewGrid (new_grid, NetCdfFileName, stdName='air_temperature', time = ncTime) !``` ! ! Read grid from _ESRI ASCII_ file: !```fortran ! TYPE(grid_real) :: new_grid ! fileName = 'file.asc' ! CALL NewGrid (new_grid, fileName, ESRI_ASCII) !``` ! ! Read grid from _ESRI BINARY_ file: !```fortran ! TYPE(grid_real) :: new_grid ! fileName = 'file' ! CALL NewGrid (new_grid, fileName, ESRI_BINARY) !``` ! ! Initialize grid using an existing one as template. The new grid ! is initialized with value 0: !```fortran ! TYPE(grid_real) :: template ! TYPE(grid_integer) :: new_Grid ! CALL NewGrid (newGrid, template) !``` ! If you want assign an initial value different than 0: !```fortran ! CALL NewGrid (new_Grid, template, 100) !``` ! The initial value must be of the same type as values stored in grid ! ! To export a grid to a file you can: ! ! write the grid to a netcdf file as variable termed VariableName: !```fortran ! TYPE(grid_real) :: grid ! ... some statements that assign data to grid ! CALL ExportGrid (grid, NetCdfFileName, VariableName, 'add') !``` ! If the netcdf file is multitemporal, you can add grid to already existing ones in file: !```fortran ! TYPE(grid_real) :: grid ! ... some statements that assign data to grid ! CALL ExportGrid (grid, NetCdfFileName, VariableName, 'append') !``` ! If NetCdfFileName does not exist, it is created. ! ! Write the grid to _ESRI ASCII_ or _BINARY_ file: !```fortran ! TYPE(grid_real) :: grid ! ... some statements that assign data to grid ! fileName = 'file.asc' ! CALL ExportGrid (grid, fileName, ESRI_ASCII) ! fileName = 'file' ! CALL ExportGrid (grid, fileName, ESRI_BINARY) !``` ! If a grid is no more used you can free memory: !```fortran ! TYPE(grid_real) :: grid ! ... some statements that assign data to grid ! CALL GridDestroy (grid) !``` ! ! References: ! NetCDF: http://www.unidata.ucar.edu/software/netcdf/ ! CF 1.0: http://cf-pcmdi.llnl.gov/ ! UDUNITS: http://www.unidata.ucar.edu/software/udunits/ MODULE GridLib ! Modules used: USE netcdf , ONLY : & ! Imported Parameters: NF90_NOWRITE, nf90_noerr, NF90_CLOBBER, & NF90_FLOAT, NF90_DOUBLE, NF90_GLOBAL, NF90_WRITE, & NF90_UNLIMITED, NF90_INT, NF90_MAX_VAR_DIMS, & ! Imported Routines: nf90_open, nf90_close, nf90_inq_dimid, nf90_inq_varid, & nf90_inquire_dimension, nf90_get_var, nf90_strerror, & nf90_inquire, nf90_get_att, nf90_inquire_variable, & nf90_create, nf90_def_dim, nf90_def_var, nf90_put_att, & nf90_enddef, nf90_put_var, nf90_redef USE DataTypeSizes ,ONLY: & ! Imported Parameters: short,long,float,double USE LogLib, ONLY : & ! Imported Routines: Catch USE ErrorCodes, ONLY : & ! Imported parameters: ncIOError, memAllocError, openFileError, genericIOError, & unknownOption USE Chronos, ONLY : & ! Imported type definitions: DateTime , & ! Imported routines: UtcNow, ToUtc, DateTimeIsDefault, & ASSIGNMENT(=), OPERATOR(-), & OPERATOR(+), OPERATOR(==), & OPERATOR(>), OPERATOR(<), & !Imported parameters: timeDefault, timeString USE GeoLib, ONLY : & ! Imported type definitions: CRS, & ! Imported parameters: GEODETIC, UTM, GAUSS_BOAGA, TM, & WGS84, ED50, ROME40, & ! Imported routines: SetCRS, SetTransverseMercatorParameters, SetGeodeticParameters IMPLICIT NONE ! Global (i.e. public) Declarations: ! Global Type Definitions: TYPE grid_real INTEGER (KIND = short) :: jdim !!number of columns INTEGER (KIND = short) :: idim !!number of rows CHARACTER (LEN = 300) :: standard_name !! CF 1.0 compliant standard name CHARACTER (LEN = 300) :: var_name !! name of the variable CHARACTER (LEN = 300) :: long_name !! long descriptive name CHARACTER (LEN = 300) :: file_name !! name of the file from which grid was read CHARACTER (LEN = 30) :: units !!UDUNITS compliant measure units CHARACTER (LEN = 20) :: varying_mode !!mode to vary: sequence, linear REAL (KIND = float) :: nodata !!scalar identifying missing value REAL (KIND = float) :: valid_min !!minimum valid value REAL (KIND = float) :: valid_max !!maximum valid value TYPE (DateTime) :: reference_time !!reference time from which calculate current TYPE (DateTime) :: current_time !!current date and time of the grid in memory TYPE (DateTime) :: next_time !!time when next update is required INTEGER (KIND = short) :: time_index !!position of grid in time dimension in netcdf file CHARACTER (LEN = 7) :: time_unit !!define time unit. Accepted values are: !!seconds, second, sec, s !!minutes, minute, min !!hours, hour, hr, h !!days, day, d REAL (KIND = float), ALLOCATABLE :: mat(:,:) !!data variable contained in grid INTEGER (KIND = short), ALLOCATABLE :: Iraster (:,:) !!used to transform not regular grid to regular raster INTEGER (KIND = short), ALLOCATABLE :: Jraster (:,:) !!used to transform not regular grid to regular raster !!georeferencing informations REAL (KIND = float) :: cellsize REAL (KIND = float) :: xllcorner REAL (KIND = float) :: yllcorner CHARACTER (LEN = 1000) :: esri_pe_string !!used by ArcMap 9.2 TYPE (CRS) :: grid_mapping !!Coordinate reference System END TYPE grid_real TYPE grid_integer INTEGER (KIND = short) :: jdim !!number of columns INTEGER (KIND = short) :: idim !!number of rows CHARACTER (LEN = 300) :: standard_name !! CF 1.0 compliant standard name CHARACTER (LEN = 300) :: var_name !! name of the variable CHARACTER (LEN = 300) :: long_name !! long descriptive name CHARACTER (LEN = 300) :: file_name !! name of the file from which grid was read CHARACTER (LEN = 30) :: units !!UDUNITS compliant measure units CHARACTER (LEN = 20) :: varying_mode !!mode to vary: steady, sequence, linear INTEGER (KIND = short) :: nodata !!scalar identifying missing value INTEGER (KIND = long) :: valid_min !!minimum valid value INTEGER (KIND = long) :: valid_max !!maximum valid value TYPE (DateTime) :: reference_time !!reference time from which calculate current TYPE (DateTime) :: current_time !!current date and time of the grid in memory TYPE (DateTime) :: next_time !!time when next update is required INTEGER (KIND = short) :: time_index !!position of grid in time dimension in netcdf file CHARACTER (LEN = 7) :: time_unit !!define time unit. Accepted values are: !!seconds, second, sec, s !!minutes, minute, min !!hours, hour, hr, h !!days, day, d INTEGER (KIND = long), ALLOCATABLE :: mat(:,:)!!data contained in grid INTEGER (KIND = short), ALLOCATABLE :: Iraster (:,:) !!used to transform not regular grid to regular raster INTEGER (KIND = short), ALLOCATABLE :: Jraster (:,:) !!used to transform not regular grid to regular raster !!georeferencing informations REAL (KIND = float) :: cellsize REAL (KIND = float) :: xllcorner REAL (KIND = float) :: yllcorner CHARACTER (LEN = 1000) :: esri_pe_string !!used by ArcMap 9.2 TYPE (CRS) :: grid_mapping !!Coordinate reference System END TYPE grid_integer ! Global Parameters: INTEGER (KIND = short), PARAMETER :: ESRI_ASCII = 1 INTEGER (KIND = short), PARAMETER :: ESRI_BINARY = 2 INTEGER (KIND = short), PARAMETER :: NET_CDF = 3 ! Global Scalars: ! Global Arrays: !Global Procedures: PUBLIC :: Newgrid PUBLIC :: ExportGrid PUBLIC :: GridDestroy PUBLIC :: SetStandardName PUBLIC :: SetLongName PUBLIC :: SetUnits PUBLIC :: SetVaryingMode PUBLIC :: SetReferenceTime PUBLIC :: SetCurrentTime PUBLIC :: SetEsriPeString PUBLIC :: GetGridMapping PUBLIC :: GetDtGrid PUBLIC :: GetFirstDate PUBLIC :: GetTimeSteps ! Local (i.e. private) Declarations: ! Local Procedures: PRIVATE :: ParseTime PRIVATE :: TimeIndex PRIVATE :: NewGridFloatFromNetCDF PRIVATE :: NewGridIntegerFromNetCDF PRIVATE :: NewGridFloatFromFile PRIVATE :: NewGridFloatFromESRI_ASCII PRIVATE :: NewGridFloatFromESRI_BINARY PRIVATE :: NewGridIntegerFromFile PRIVATE :: NewGridIntegerFromESRI_ASCII PRIVATE :: NewGridIntegerFromESRI_BINARY PRIVATE :: NewGridFloatAsGridFloat PRIVATE :: NewGridFloatAsGridInteger PRIVATE :: NewGridIntegerAsGridInteger PRIVATE :: NewGridIntegerAsGridFloat PRIVATE :: GridDestroyFloat PRIVATE :: GridDestroyInteger PRIVATE :: GetXYSizesFromFile PRIVATE :: GetGeoreferenceFromNCdataSet PRIVATE :: ExportGridFloatToNetCDF PRIVATE :: ExportGridFloatToFile PRIVATE :: ExportGridFloatToESRI_ASCII PRIVATE :: ExportGridFloatToESRI_BINARY PRIVATE :: ExportGridIntegerToFile PRIVATE :: ExportGridIntegerToESRI_ASCII PRIVATE :: ExportGridIntegerToESRI_BINARY PRIVATE :: SetStandardNameFloat PRIVATE :: SetStandardNameInteger PRIVATE :: SetLongNameFloat PRIVATE :: SetLongNameInteger PRIVATE :: SetUnitsFloat PRIVATE :: SetUnitsInteger PRIVATE :: SetVaryingModeFloat PRIVATE :: SetVaryingModeInteger PRIVATE :: SetReferenceTimeFloat PRIVATE :: SetReferenceTimeInteger PRIVATE :: SetCurrentTimeFloat PRIVATE :: SetCurrentTimeInteger PRIVATE :: SetEsriPeStringFloat PRIVATE :: SetEsriPeStringInteger PRIVATE :: SwapGridRealForward PRIVATE :: SwapGridIntegerForward PRIVATE :: SwapGridRealBack PRIVATE :: SwapGridIntegerBack PRIVATE :: NextTime PRIVATE :: GetTimeStepsFromFile PRIVATE :: GetTimeStepsFromNCid PRIVATE :: ScanTimeStringAsString ! Local Type Definitions: ! Local Parameters: INTEGER (KIND = short), PRIVATE, PARAMETER :: MISSING_DEF_INT = -9999 REAL (KIND = float), PRIVATE, PARAMETER :: MISSING_DEF_REAL = -9999.9 ! Local Scalars: ! Local Arrays: ! Operator definitions: ! Define new operators or overload existing ones. INTERFACE NewGrid !MODULE PROCEDURE NewGridFloatFromNetCDF !MODULE PROCEDURE NewGridIntegerFromNetCDF MODULE PROCEDURE NewGridFloatFromFile MODULE PROCEDURE NewGridIntegerFromFile MODULE PROCEDURE NewGridFloatAsGridFloat MODULE PROCEDURE NewGridFloatAsGridInteger MODULE PROCEDURE NewGridIntegerAsGridInteger MODULE PROCEDURE NewGridIntegerAsGridFloat END INTERFACE INTERFACE ExportGrid MODULE PROCEDURE ExportGridFloatToNetCDF MODULE PROCEDURE ExportGridIntegerToNetCDF MODULE PROCEDURE ExportGridFloatToFile MODULE PROCEDURE ExportGridIntegerToFile END INTERFACE INTERFACE GridDestroy MODULE PROCEDURE GridDestroyFloat MODULE PROCEDURE GridDestroyInteger END INTERFACE INTERFACE GetTimeSteps MODULE PROCEDURE GetTimeStepsFromNCid MODULE PROCEDURE GetTimeStepsFromFile END INTERFACE INTERFACE SetStandardName MODULE PROCEDURE SetStandardNameFloat MODULE PROCEDURE SetStandardNameInteger END INTERFACE INTERFACE SetLongName MODULE PROCEDURE SetLongNameFloat MODULE PROCEDURE SetLongNameInteger END INTERFACE INTERFACE SetUnits MODULE PROCEDURE SetUnitsFloat MODULE PROCEDURE SetUnitsInteger END INTERFACE INTERFACE SetVaryingMode MODULE PROCEDURE SetVaryingModeFloat MODULE PROCEDURE SetVaryingModeInteger END INTERFACE INTERFACE SetReferenceTime MODULE PROCEDURE SetReferenceTimeFloat MODULE PROCEDURE SetReferenceTimeInteger END INTERFACE INTERFACE SetCurrentTime MODULE PROCEDURE SetCurrentTimeFloat MODULE PROCEDURE SetCurrentTimeInteger END INTERFACE INTERFACE SetEsriPeString MODULE PROCEDURE SetEsriPeStringFloat MODULE PROCEDURE SetEsriPeStringInteger END INTERFACE INTERFACE GetGridMapping MODULE PROCEDURE GetGridMappingFloat MODULE PROCEDURE GetGridMappingInteger END INTERFACE !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! create a new raster_real reading data from NetCDF file ! The variable to read can be defined by its current name or ! the standard_name. The dimensions x and y of the variable ! is calculated searching from the dimensions of the couple of ! variables with 'standard_name' equal to 'projection_x_coordinate' ! and 'projection_y_coordinate' for projected reference systems ! or 'longitude' and 'latitude' for geographic reference systems ! or 'grid_longitude' and 'grid_latitude' for rotated pole systems ! If a comprehensible reference systems is not found a geodetic ! reference system is supposed. ! Once the variable is retrieved, offset and scale factor are applied ! and a check on minimum and maximum valid value is performed. SUBROUTINE NewGridFloatFromNetCDF & ! (layer, fileName, variable, stdName, time) USE StringManipulation, ONLY: & !Imported routines: StringCompact IMPLICIT NONE ! Arguments with intent(in): CHARACTER (LEN = *), INTENT(in) :: fileName !!NetCDF file to be read CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: variable !!variable to read CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: stdName !!standard name of !!the variable to read TYPE (DateTime), OPTIONAL, INTENT(in) :: time !!time of the grid to read ! Arguments with intent(out): TYPE (grid_real), INTENT (out) :: layer !!gridreal to return ! Local scalars: INTEGER (KIND = short) :: ios !!error return code INTEGER (KIND = short) :: ncStatus !!error code returned by NetCDF routines INTEGER (KIND = short) :: ncId !!NetCdf Id for the file INTEGER (KIND = short) :: varId !!variable Id INTEGER (KIND = short) :: nVars !!number of variables CHARACTER (LEN = 80) :: attribute INTEGER (KIND = short) :: nDimsVar !!number of dimensions of a variable INTEGER (KIND = short) :: i, j !!loop index CHARACTER (LEN = 100) :: variableName CHARACTER (LEN = 1) :: shp REAL (KIND = float) :: scale_factor REAL (KIND = float) :: offset INTEGER (KIND = short) :: latlon ! Local arrays: INTEGER, ALLOCATABLE :: slice (:) REAL (KIND = float), ALLOCATABLE :: tempGrid (:,:) !------------end of declaration------------------------------------------------ !------------------------------------------------------------------------------ ![1.0] open NetCDF dataset with read-only access !------------------------------------------------------------------------------ ncStatus = nf90_open (fileName, NF90_NOWRITE, ncId) IF (ncStatus /= nf90_noerr) THEN CALL Catch ('error', 'GridLib', & TRIM (nf90_strerror (ncStatus) )//': ', & code = ncIOError, argument = fileName ) ENDIF !------------------------------------------------------------------------------ ![2.0] get x and y size and allocate array !------------------------------------------------------------------------------ CALL GetXYSizesFromFile (ncId, layer % jdim, layer % idim, shape = shp, latlon = latlon) !allocate grid ALLOCATE ( layer % mat (layer % idim, layer % jdim), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError,argument = variable ) ENDIF !allocate temporary grid IF (latlon == 1) THEN !coordinate order lon-lat ALLOCATE ( tempGrid (layer % jdim, layer % idim), STAT = ios ) ELSE !coordinate order lat-lon ALLOCATE ( tempGrid (layer % idim, layer % jdim), STAT = ios ) END IF IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError,argument = variable ) ENDIF !------------------------------------------------------------------------------ ![3.0] get values !------------------------------------------------------------------------------ IF ( PRESENT (variable) ) THEN ncStatus = nf90_inq_varid (ncId, variable, varId) CALL ncErrorHandler (ncStatus) ELSE IF (PRESENT (stdName) ) THEN !search variable corresponding to standard name !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nVariables = nVars ) CALL ncErrorHandler (ncStatus) DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN !'standard_name' is found IF ( StringCompact(attribute(1:LEN_TRIM(attribute))) == stdName(1:LEN_TRIM(stdName)) )THEN varId = i END IF END IF END DO END IF !verify number of dimensions of variable to read ncstatus = nf90_inquire_variable(ncId, varId, ndims = nDimsVar) CALL ncErrorHandler (ncStatus) !If number of dimensions = 3 read info about time IF (nDimsVar == 3) THEN ALLOCATE ( slice(3) ) !Read time information CALL ParseTime (ncId, layer % reference_time, layer % time_unit) !Define current time slice (if present) IF (PRESENT (time) ) THEN slice = (/ 1, 1, 1 /) slice(3) = TimeIndex (ncId, layer % reference_time, & layer % time_unit, time) layer % current_time = time layer % time_index = slice(3) ELSE !read the first grid slice = (/ 1, 1, 1 /) layer % time_index = 1 !set current time layer % current_time = TimeByIndex (ncId, layer % reference_time, & layer % time_unit, layer % time_index) ENDIF !define time of the next grid layer % next_time = NextTime (ncId, layer % reference_time, & layer % time_unit, layer % time_index) ELSE ALLOCATE ( slice(2) ) slice = (/ 1, 1/) layer % current_time = timeDefault layer % next_time = timeDefault END IF ncStatus = nf90_get_var (ncId, varId, tempGrid , start = slice) CALL ncErrorHandler (ncStatus) !transpose temporary matrix to grid specification CALL SwapGridRealForward (tempGrid, layer % mat, latlon) !deallocate temporary grid DEALLOCATE (tempGrid) !------------------------------------------------------------------------------ ![4.0] get attributes !------------------------------------------------------------------------------ ncStatus = nf90_get_att (ncId, varId, name = 'standard_name', & values = layer % standard_name) ncStatus = nf90_get_att (ncId, varId, name = 'long_name', & values = layer % long_name) ncStatus = nf90_get_att (ncId, varId, name = 'units', & values = layer % units) ncStatus = nf90_get_att (ncId, varId, name = 'varying_mode', & values = layer % varying_mode) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) layer % varying_mode = 'sequence' ncStatus = nf90_get_att (ncId, varId, name = '_FillValue', & values = layer % nodata) !if _FillValue is not defined search for missing_value IF (ncStatus /= nf90_noerr) THEN ncStatus = nf90_get_att (ncId, varId, name = 'missing_value', & values = layer % nodata) END IF !if missing_value is not defined set to default IF (ncStatus /= nf90_noerr) THEN layer % nodata = MISSING_DEF_REAL END IF ncStatus = nf90_get_att (ncId, varId, name = 'valid_min', & values = layer % valid_min) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) layer % valid_min = layer % nodata ncStatus = nf90_get_att (ncId, varId, name = 'valid_max', & values = layer % valid_max) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) layer % valid_max = layer % nodata ncStatus = nf90_get_att (ncId, varId, name = 'scale_factor', & values = scale_factor) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) scale_factor = 1.0 ncStatus = nf90_get_att (ncId, varId, name = 'add_offset', & values = offset) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) offset = 0. ncStatus = nf90_get_att (ncId, varId, name = 'esri_pe_string', & values = layer % esri_pe_string) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) layer % esri_pe_string = '' !------------------------------------------------------------------------------ ![5.0] check values and apply scale factor and offset !------------------------------------------------------------------------------ !retrieve name of variable ncstatus = nf90_inquire_variable(ncId, varId, name = variableName) layer % var_name = variableName DO i = 1, layer % idim DO j = 1, layer % jdim IF ( layer % mat (i,j) /= layer % nodata ) THEN !apply scale factor layer % mat (i,j) = layer % mat (i,j) * scale_factor !add offset layer % mat (i,j) = layer % mat (i,j) + offset !check lower bound IF ( layer % valid_min /= layer % nodata .AND. & layer % mat (i,j) < layer % valid_min ) THEN layer % mat (i,j) = layer % valid_min CALL Catch ('info', 'GridLib', 'corrected value exceeding & lower bound in variable: ', argument = variableName ) END IF !check upper bound IF ( layer % valid_max /= layer % nodata .AND. & layer % mat (i,j) > layer % valid_max ) THEN layer % mat (i,j) = layer % valid_max CALL Catch ('info', 'GridLib', 'corrected value exceeding & upper bound in variable: ', argument = variableName ) END IF END IF END DO END DO !------------------------------------------------------------------------------ ![6.0] set file name !------------------------------------------------------------------------------ layer % file_name = fileName !------------------------------------------------------------------------------ ![7.0] read georeferencing informations from netCDF file ! Regularize grid if necessary !------------------------------------------------------------------------------ IF (shp == 'm') THEN !2 dimensional (matrix) coordinate for not regular grid IF (.NOT. ALLOCATED (layer % Iraster)) THEN ALLOCATE (layer % Iraster(layer%idim,layer%jdim)) ALLOCATE (layer % Jraster(layer%idim,layer%jdim)) CALL GetGeoreferenceFromNCdataSet (ncId, varId, layer % cellsize, & layer % xllcorner, layer % yllcorner, & layer % grid_mapping, layer % Iraster, layer % Jraster, .TRUE.) END IF !create a temporary copy of layer % mat IF (ALLOCATED (tempGrid) ) THEN DEALLOCATE (tempGrid) ALLOCATE (tempGrid(layer%idim,layer%jdim)) ELSE ALLOCATE (tempGrid(layer%idim,layer%jdim)) END IF tempGrid = layer % mat layer % mat = layer % nodata !fill in the regular grid DO i = 1, layer%idim DO j = 1, layer%jdim IF (layer % Iraster (i,j) /= -9999 .AND. layer % Jraster (i,j) /= -9999) THEN layer % mat (layer%Iraster(i,j),layer%Jraster(i,j)) = tempGrid (i,j) ELSE layer % mat(i,j) = layer % nodata END IF END DO END DO DEALLOCATE (tempGrid) ELSE !regular grid stores coordinate in 1dimensional array (vector) CALL GetGeoreferenceFromNCdataSet (ncId, varId, layer % cellsize, & layer % xllcorner, layer % yllcorner, & layer % grid_mapping, layer % Iraster, layer % Iraster, .FALSE.) END IF !------------------------------------------------------------------------------ ![8.0] close NetCDF dataset !------------------------------------------------------------------------------ ncStatus = nf90_close (ncid) CALL ncErrorHandler (ncStatus) END SUBROUTINE NewGridFloatFromNetCDF !============================================================================== !| Description: ! create a new grid_integer reading data from NetCDF file ! The variable to read can be defined by its current name or ! the standard_name. The dimensions x and j of the variable ! is calculated searching from the dimensions of the couple of ! variables with 'standard_name' equal to 'projection_x_coordinate' ! and 'projection_y_coordinate' for projected reference systems ! or 'longitude' and 'latitude' for geographic reference systems ! or 'grid_longitude' and 'grid_latitude' for rotated pole systems ! If a comprehensible reference systems is not found a geodetic ! reference system is supposed. ! Once the variable is retrieved, offset and scale factor are applied ! and a check on minimum and maximum valid value is performed. SUBROUTINE NewGridIntegerFromNetCDF & ! (layer, fileName, variable, stdName, time) USE StringManipulation, ONLY: & !Imported routines: StringCompact IMPLICIT NONE ! Scalar arguments with intent(in): CHARACTER (LEN = *), INTENT(in) :: fileName !!NetCDF file to read CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: variable !!variable to read CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: stdName !!standard name of !!the variable to read TYPE (DateTime), OPTIONAL, INTENT(in) :: time !!time of the grid to read ! Arguments with intent(out): TYPE (grid_integer), INTENT (out) :: layer !!gridreal to return ! Local scalars: INTEGER (KIND = short) :: ios !!error return code INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: ncId !!NetCdf Id for the file INTEGER (KIND = short) :: varId !!variable Id INTEGER (KIND = short) :: nVars !!number of variables CHARACTER (LEN = 80) :: attribute INTEGER (KIND = short) :: nDimsVar !!number of dimensions of a variable INTEGER (KIND = short) :: i, j !!loop index CHARACTER (LEN = 100) :: variableName CHARACTER (LEN = 1) :: shp INTEGER (KIND = long) :: scale_factor INTEGER (KIND = long) :: offset ! Local arrays: INTEGER, ALLOCATABLE :: slice (:) INTEGER (KIND = long), ALLOCATABLE :: tempGrid (:,:) !------------end of declaration------------------------------------------------ !------------------------------------------------------------------------------ ![1.0] open NetCDF dataset with read-only access !------------------------------------------------------------------------------ ncStatus = nf90_open (fileName, NF90_NOWRITE, ncId) IF (ncStatus /= nf90_noerr) THEN CALL Catch ('error', 'GridLib', & TRIM (nf90_strerror (ncStatus) )//': ', & code = ncIOError, argument = fileName ) ENDIF !------------------------------------------------------------------------------ ![2.0] get x and y size and allocate array !------------------------------------------------------------------------------ CALL GetXYSizesFromFile (ncId, layer % jdim, layer % idim, shape = shp) !allocate grid ALLOCATE ( layer % mat (layer % idim, layer % jdim), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError,argument = variable ) ENDIF !allocate temporary grid ALLOCATE ( tempGrid (layer % jdim, layer % idim), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError,argument = variable ) ENDIF !------------------------------------------------------------------------------ ![3.0] get values !------------------------------------------------------------------------------ IF ( PRESENT (variable) ) THEN ncStatus = nf90_inq_varid (ncId, variable, varId) CALL ncErrorHandler (ncStatus) ELSE IF (PRESENT (stdName) ) THEN !search variable corresponding to standard name !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nVariables = nVars ) CALL ncErrorHandler (ncStatus) DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN !error is raised if 'standard_name' is not found IF ( StringCompact(attribute(1:LEN_TRIM(attribute))) == stdName(1:LEN_TRIM(stdName)) )THEN varId = i END IF END IF END DO END IF !verify number of dimensions of variable to read ncstatus = nf90_inquire_variable(ncId, varId, ndims = nDimsVar) CALL ncErrorHandler (ncStatus) !If number of dimensions = 3 read info about time IF (nDimsVar == 3) THEN ALLOCATE ( slice(3) ) !Read time informations CALL ParseTime (ncId, layer % reference_time, layer % time_unit) !Define current time slice (if present) IF (PRESENT (time) ) THEN slice = (/ 1, 1, 1 /) slice(3) = TimeIndex (ncId, layer % reference_time, & layer % time_unit, time) layer % current_time = time layer % time_index = slice(3) ELSE slice = (/ 1, 1, 1 /) ENDIF ELSE ALLOCATE ( slice(2) ) slice = (/ 1, 1/) END IF ncStatus = nf90_get_var (ncId, varId, tempGrid , start = slice) CALL ncErrorHandler (ncStatus) !transpose temporary matrix to grid specification CALL SwapGridIntegerForward (tempGrid, layer % mat) !deallocate temporary grid DEALLOCATE (tempGrid) !------------------------------------------------------------------------------ ![4.0] get attributes !------------------------------------------------------------------------------ ncStatus = nf90_get_att (ncId, varId, name = 'standard_name', & values = layer % standard_name) ncStatus = nf90_get_att (ncId, varId, name = 'long_name', & values = layer % long_name) ncStatus = nf90_get_att (ncId, varId, name = 'units', & values = layer % units) ncStatus = nf90_get_att (ncId, varId, name = 'varying_mode', & values = layer % varying_mode) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) layer % varying_mode = 'sequence' ncStatus = nf90_get_att (ncId, varId, name = '_FillValue', & values = layer % nodata) !if _FillValue is not defined search for missing_value IF (ncStatus /= nf90_noerr) THEN ncStatus = nf90_get_att (ncId, varId, name = 'missing_value', & values = layer % nodata) END IF !if missing_value is not defined set to default IF (ncStatus /= nf90_noerr) THEN layer % nodata = MISSING_DEF_REAL END IF ncStatus = nf90_get_att (ncId, varId, name = 'valid_min', & values = layer % valid_min) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) layer % valid_min = layer % nodata ncStatus = nf90_get_att (ncId, varId, name = 'valid_max', & values = layer % valid_max) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) layer % valid_max = layer % nodata ncStatus = nf90_get_att (ncId, varId, name = 'scale_factor', & values = scale_factor) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) scale_factor = 1 ncStatus = nf90_get_att (ncId, varId, name = 'add_offset', & values = offset) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) offset = 0 ncStatus = nf90_get_att (ncId, varId, name = 'esri_pe_string', & values = layer % esri_pe_string) !If attribute is not defined set to default IF (ncStatus /= nf90_noerr) layer % esri_pe_string = '' !------------------------------------------------------------------------------ ![5.0] check values and apply scale factor and offset !------------------------------------------------------------------------------ !retrieve name of variable ncstatus = nf90_inquire_variable(ncId, varId, name = variableName) layer % var_name = variableName DO i = 1, layer % jdim DO j = 1, layer % idim IF ( layer % mat (i,j) /= layer % nodata ) THEN !apply scale factor layer % mat (i,j) = layer % mat (i,j) * scale_factor !add offset layer % mat (i,j) = layer % mat (i,j) + offset !check lower bound IF ( layer % valid_min /= layer % nodata .AND. & layer % mat (i,j) < layer % valid_min ) THEN layer % mat (i,j) = layer % valid_min CALL Catch ('info', 'GridLib', 'corrected value exceeding & lower bound in variable: ', argument = variableName ) END IF !check upper bound IF ( layer % valid_max /= layer % nodata .AND. & layer % mat (i,j) > layer % valid_max ) THEN layer % mat (i,j) = layer % valid_max CALL Catch ('info', 'GridLib', 'corrected value exceeding & upper bound in variable: ', argument = variableName ) END IF END IF END DO END DO !------------------------------------------------------------------------------ ![6.0] set file name !------------------------------------------------------------------------------ layer % file_name = fileName !------------------------------------------------------------------------------ ![7.0] read georeferencing informations from netCDF file !------------------------------------------------------------------------------ IF (shp == 'm') THEN !2 dimensional (matrix) coordinate for not regular grid IF (.NOT. ALLOCATED (layer % Iraster)) THEN ALLOCATE (layer % Iraster(layer%idim,layer%jdim)) ALLOCATE (layer % Jraster(layer%idim,layer%jdim)) CALL GetGeoreferenceFromNCdataSet (ncId, varId, layer % cellsize, & layer % xllcorner, layer % yllcorner, & layer % grid_mapping, layer % Iraster, layer % Jraster, .TRUE.) END IF !create a temporary copy of layer % mat IF (ALLOCATED (tempGrid) ) THEN DEALLOCATE (tempGrid) ALLOCATE (tempGrid(layer%idim,layer%jdim)) ELSE ALLOCATE (tempGrid(layer%idim,layer%jdim)) END IF tempGrid = layer % mat layer % mat = layer % nodata !fill in the regular grid DO i = 1, layer%idim DO j = 1, layer%jdim IF (layer % Iraster (i,j) /= -9999 .AND. layer % Jraster (i,j) /= -9999) THEN layer % mat (layer%Iraster(i,j),layer%Jraster(i,j)) = tempGrid (i,j) ELSE layer % mat(i,j) = layer % nodata END IF END DO END DO DEALLOCATE (tempGrid) ELSE !regular grid stores coordinate in 1dimensional array (vector) CALL GetGeoreferenceFromNCdataSet (ncId, varId, layer % cellsize, & layer % xllcorner, layer % yllcorner, & layer % grid_mapping, layer % Iraster, layer % Iraster, .FALSE.) END IF !------------------------------------------------------------------------------ ![8.0] close NetCDF dataset !------------------------------------------------------------------------------ ncStatus = nf90_close (ncid) CALL ncErrorHandler (ncStatus) END SUBROUTINE NewGridIntegerFromNetCDF !============================================================================== !| Description: ! read a grid from a file. ! List of supported format: ! ESRI_ASCII: ESRI ASCII GRID ! ESRI_BINARY: ESRI BINARY GRID ! NET_CDF: NetCDF CF compliant SUBROUTINE NewGridFloatFromFile & ! (layer, fileName, fileFormat, variable, stdName, time) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: fileName !! file to read INTEGER (KIND = short), INTENT(IN) :: fileFormat !! format of the file to read CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: variable !!variable to read CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: stdName !!standard name of !!the variable to read TYPE (DateTime), OPTIONAL, INTENT(in) :: time !!time of the grid to read !Arguments with intent(out): TYPE (grid_real), INTENT(OUT) :: layer !!gridreal to return !Local variables: !------------end of declaration------------------------------------------------ IF ( fileformat == ESRI_ASCII ) THEN CALL NewGridFloatFromESRI_ASCII (fileName, layer) ELSE IF ( fileformat == ESRI_BINARY ) THEN CALL NewGridFloatFromESRI_BINARY (fileName, layer) ELSE IF ( fileformat == NET_CDF ) THEN IF (PRESENT(stdName)) THEN IF (PRESENT (time)) THEN CALL NewGridFloatFromNetCDF (layer, fileName, stdName = stdName, time = time) ELSE CALL NewGridFloatFromNetCDF (layer, fileName, stdName= stdName) END IF ELSE IF (PRESENT(variable)) THEN IF (PRESENT (time)) THEN CALL NewGridFloatFromNetCDF (layer, fileName, variable = variable, time = time) ELSE CALL NewGridFloatFromNetCDF (layer, fileName, variable = variable) END IF END IF ELSE CALL Catch ('error', 'GridLib', & 'unknown option in reading file grid: ', & code = unknownOption, argument = fileName ) END IF END SUBROUTINE NewGridFloatFromFile !============================================================================== !! Description: !! read a float grid from a ESRI ASCII file. SUBROUTINE NewGridFloatFromESRI_ASCII & ! (fileName, layer) USE Utilities, ONLY: & !Imported routines: GetUnit IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(in) :: fileName !! file to read !Arguments with intent(out) TYPE (grid_real), INTENT (OUT) :: layer !!returned grid float !Local variables: INTEGER (KIND = short) :: fileUnit INTEGER (KIND = short) :: ios CHARACTER (LEN = 50) :: dummy INTEGER (KIND = short) :: i, j REAL (KIND = double) :: corner !------------end of declaration------------------------------------------------ !open file fileUnit = GetUnit () OPEN (UNIT = fileUnit, file = fileName, STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = fileName ) END IF !read number of columns READ (fileUnit,*,IOSTAT = ios) dummy, layer % jdim IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading columns in file: ', & code = genericIOError, argument = fileName ) END IF !read number of rows READ (fileUnit,*,IOSTAT = ios) dummy, layer % idim IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading rows in file: ', & code = genericIOError, argument = fileName ) END IF !read xll corner READ (fileUnit,*,IOSTAT = ios) dummy, corner layer % xllcorner = corner IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading xll in file: ', & code = genericIOError, argument = fileName ) END IF !read yll corner READ (fileUnit,*,IOSTAT = ios) dummy, corner layer % yllcorner = corner IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading yll in file: ', & code = genericIOError, argument = fileName ) END IF !read cellsize READ (fileUnit,*,IOSTAT = ios) dummy, layer % cellsize IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading cellsize in file: ', & code = genericIOError, argument = fileName ) END IF !read nodata value READ (fileUnit,*,IOSTAT = ios) dummy, layer % nodata IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading missing value in file: ', & code = genericIOError, argument = fileName ) END IF !allocate grid ALLOCATE ( layer % mat (layer % idim, layer % jdim), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError,argument = fileName ) ENDIF !read data DO i = 1,layer % idim READ (fileUnit,*) ( layer % mat (i,j) , j = 1,layer % jdim ) END DO CLOSE (fileUnit) !Set to default other fields layer % standard_name = '' layer % long_name = '' layer % units = '' layer % varying_mode = 'sequence' layer % valid_min = layer % nodata layer % valid_max = layer % nodata layer % esri_pe_string = '' layer % reference_time = timeDefault layer % current_time = timeDefault layer % next_time = timeDefault END SUBROUTINE NewGridFloatFromESRI_ASCII !============================================================================== !| Description: ! read a float grid from a ESRI BINARY file. SUBROUTINE NewGridFloatFromESRI_BINARY & ! (fileName, layer) USE Utilities, ONLY: & !Imported routines: GetUnit IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(in) :: fileName !! file to read !Arguments with intent(out) TYPE (grid_real), INTENT (OUT) :: layer !!returned grid float !Local variables: INTEGER (KIND = short) :: fileUnit INTEGER (KIND = short) :: ios CHARACTER (LEN = 50) :: dummy INTEGER (KIND = short) :: recordLength INTEGER (KIND = short) :: recordNumber INTEGER (KIND = short) :: i, j !------------end of declaration------------------------------------------------ !open file fileUnit = GetUnit () OPEN (UNIT = fileUnit, file = TRIM(fileName) // '.hdr', & STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = TRIM(fileName) // '.hdr' ) END IF !read number of columns READ (fileUnit,*,IOSTAT = ios) dummy, layer % jdim IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading columns in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read number of rows READ (fileUnit,*,IOSTAT = ios) dummy, layer % idim IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading rows in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read xll corner READ (fileUnit,*,IOSTAT = ios) dummy, layer % xllcorner IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading xll in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read yll corner READ (fileUnit,*,IOSTAT = ios) dummy, layer % yllcorner IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading yll in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read cellsize READ (fileUnit,*,IOSTAT = ios) dummy, layer % cellsize IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading cellsize in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read nodata value READ (fileUnit,*,IOSTAT = ios) dummy, layer % nodata IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading missing value in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !allocate grid ALLOCATE ( layer % mat (layer % idim, layer % jdim), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError,argument = fileName ) ENDIF CLOSE (fileUnit) fileUnit = GetUnit () INQUIRE (IOLENGTH = recordLength) 100_float OPEN (UNIT = fileUnit, file = TRIM(fileName) // '.flt', & FORM = 'UNFORMATTED', ACCESS = 'DIRECT', RECL = recordLength, & STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = TRIM(fileName) // '.flt' ) END IF !read data recordNumber = 0 DO i = 1,layer % idim DO j = 1, layer % jdim recordNumber = recordNumber + 1 READ (fileUnit,REC = recordNumber) layer % mat (i,j) END DO END DO CLOSE (fileUnit) !Set other fields to default layer % standard_name = '' layer % long_name = '' layer % units = '' layer % varying_mode = 'sequence' layer % valid_min = layer % nodata layer % valid_max = layer % nodata layer % esri_pe_string = '' layer % reference_time = timeDefault layer % current_time = timeDefault layer % next_time = timeDefault END SUBROUTINE NewGridFloatFromESRI_BINARY !============================================================================== !| Description: ! read a grid from a file. ! ! * List of supported format: ! * ESRI_ASCII: ESRI ASCII GRID ! * ESRI_BINARY: ESRI BINARY GRID ! SUBROUTINE NewGridIntegerFromFile & ! (layer, fileName, fileFormat, variable, stdName, time) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: fileName !! file to read INTEGER (KIND = short), INTENT(IN) :: fileFormat !! format of the file to read CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: variable !!variable to read CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: stdName !!standard name of !!the variable to read TYPE (DateTime), OPTIONAL, INTENT(in) :: time !!time of the grid to read !Arguments with intent(out): TYPE (grid_integer), INTENT(OUT) :: layer !!grid to be returned !------------end of declaration------------------------------------------------ IF ( fileformat == ESRI_ASCII ) THEN CALL NewGridIntegerFromESRI_ASCII (fileName, layer) ELSE IF ( fileformat == ESRI_BINARY ) THEN CALL NewGridIntegerFromESRI_BINARY (fileName, layer) ELSE IF ( fileformat == NET_CDF ) THEN IF (PRESENT(stdName)) THEN IF (PRESENT (time)) THEN CALL NewGridIntegerFromNetCDF (layer, fileName, stdName = stdName, time = time) ELSE CALL NewGridIntegerFromNetCDF (layer, fileName, stdName= stdName) END IF ELSE IF (PRESENT(variable)) THEN IF (PRESENT (time)) THEN CALL NewGridIntegerFromNetCDF (layer, fileName, variable = variable, time = time) ELSE CALL NewGridIntegerFromNetCDF (layer, fileName, variable = variable) END IF END IF ELSE CALL Catch ('error', 'GridLib', & 'unknown option in reading file grid: ', & code = unknownOption, argument = fileName ) END IF END SUBROUTINE NewGridIntegerFromFile !============================================================================== !| Description: ! read a integer grid from a ESRI ASCII file. SUBROUTINE NewGridIntegerFromESRI_ASCII & ! (fileName, layer) USE Utilities, ONLY: & !Imported routines: GetUnit IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(in) :: fileName !! file to read !Arguments with intent(out) TYPE (grid_integer), INTENT (OUT) :: layer !!returned grid float !Local variables: INTEGER (KIND = short) :: fileUnit INTEGER (KIND = short) :: ios CHARACTER (LEN = 50) :: dummy INTEGER (KIND = short) :: i, j !------------end of declaration------------------------------------------------ !open file fileUnit = GetUnit () OPEN (UNIT = fileUnit, file = fileName, STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = fileName ) END IF !read number of columns READ (fileUnit,*,IOSTAT = ios) dummy, layer % jdim IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading columns in file: ', & code = genericIOError, argument = fileName ) END IF !read number of rows READ (fileUnit,*,IOSTAT = ios) dummy, layer % idim IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading rows in file: ', & code = genericIOError, argument = fileName ) END IF !read xll corner READ (fileUnit,*,IOSTAT = ios) dummy, layer % xllcorner IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading xll in file: ', & code = genericIOError, argument = fileName ) END IF !read yll corner READ (fileUnit,*,IOSTAT = ios) dummy, layer % yllcorner IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading yll in file: ', & code = genericIOError, argument = fileName ) END IF !read cellsize READ (fileUnit,*,IOSTAT = ios) dummy, layer % cellsize IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading cellsize in file: ', & code = genericIOError, argument = fileName ) END IF !read nodata value READ (fileUnit,*,IOSTAT = ios) dummy, layer % nodata IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading missing value in file: ', & code = genericIOError, argument = fileName ) END IF !allocate grid ALLOCATE ( layer % mat (layer % idim, layer % jdim), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError,argument = fileName ) ENDIF !read data DO i = 1,layer%idim READ (fileUnit,*) ( layer % mat (i,j) , j=1,layer % jdim ) END DO CLOSE (fileUnit) !Set to default other fields layer % standard_name = '' layer % long_name = '' layer % units = '' layer % varying_mode = 'sequence' layer % valid_min = layer % nodata layer % valid_max = layer % nodata layer % esri_pe_string = '' layer % reference_time = timeDefault layer % current_time = timeDefault layer % next_time = timeDefault END SUBROUTINE NewGridIntegerFromESRI_ASCII !============================================================================== !| Description: ! read a integer grid from a ESRI BINARY file. SUBROUTINE NewGridIntegerFromESRI_BINARY & ! (fileName, layer) USE Utilities, ONLY: & !Imported routines: GetUnit IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(in) :: fileName !! file to read !Arguments with intent(out) TYPE (grid_integer), INTENT (OUT) :: layer !!returned grid float !Local variables: INTEGER (KIND = short) :: fileUnit INTEGER (KIND = short) :: ios CHARACTER (LEN = 50) :: dummy INTEGER (KIND = short) :: recordLength INTEGER (KIND = short) :: recordNumber INTEGER (KIND = short) :: i, j !------------end of declaration------------------------------------------------ !open file fileUnit = GetUnit () OPEN (UNIT = fileUnit, file = TRIM(fileName) // '.hdr', & STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = TRIM(fileName) // '.hdr' ) END IF !read number of columns READ (fileUnit,*,IOSTAT = ios) dummy, layer % jdim IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading columns in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read number of rows READ (fileUnit,*,IOSTAT = ios) dummy, layer % idim IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading rows in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read xll corner READ (fileUnit,*,IOSTAT = ios) dummy, layer % xllcorner IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading xll in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read yll corner READ (fileUnit,*,IOSTAT = ios) dummy, layer % yllcorner IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading yll in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read cellsize READ (fileUnit,*,IOSTAT = ios) dummy, layer % cellsize IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading cellsize in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !read nodata value READ (fileUnit,*,IOSTAT = ios) dummy, layer % nodata IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error reading missing value in file: ', & code = genericIOError, argument = TRIM(fileName) // '.hdr' ) END IF !allocate grid ALLOCATE ( layer % mat (layer % idim, layer % jdim), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError,argument = fileName ) ENDIF CLOSE (fileUnit) fileUnit = GetUnit () INQUIRE (IOLENGTH = recordLength) 100_long OPEN (UNIT = fileUnit, file = TRIM(fileName) // '.flt', & FORM = 'UNFORMATTED', ACCESS = 'DIRECT', RECL = recordLength, & STATUS = 'OLD', IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = TRIM(fileName) // '.flt' ) END IF !read data recordNumber = 0 DO i = 1,layer % idim DO j = 1, layer % jdim recordNumber = recordNumber + 1 READ (fileUnit,REC = recordNumber) layer % mat (i,j) END DO END DO CLOSE (fileUnit) !Set other fields to default layer % standard_name = '' layer % long_name = '' layer % units = '' layer % varying_mode = 'sequence' layer % valid_min = layer % nodata layer % valid_max = layer % nodata layer % esri_pe_string = '' layer % reference_time = timeDefault layer % current_time = timeDefault layer % next_time = timeDefault END SUBROUTINE NewGridIntegerFromESRI_BINARY !============================================================================== !| Description: ! create a new grid_real using an existing grid_real as template SUBROUTINE NewGridFloatAsGridFloat & ! (layer, grid, initial) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_real), INTENT(in) :: grid REAL (KIND = float), OPTIONAL, INTENT(in) :: initial !Arguments with intent(out): TYPE (grid_real), INTENT(OUT) :: layer !!grid to be returned !Local variables: INTEGER (KIND = short) :: ios INTEGER (KIND = short) :: i, j !------------end of declaration------------------------------------------------ layer % jdim = grid % jdim layer % idim = grid % idim layer % varying_mode = 'sequence' !default layer % nodata = MISSING_DEF_REAL layer % valid_min = layer % nodata layer % valid_max = layer % nodata layer % cellsize = grid % cellsize layer % xllcorner = grid % xllcorner layer % yllcorner = grid % yllcorner layer % esri_pe_string = grid % esri_pe_string layer % grid_mapping = grid % grid_mapping layer % reference_time = timeDefault layer % current_time = timeDefault layer % next_time = timeDefault ALLOCATE ( layer % mat ( layer % idim, layer % jdim ), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError, argument = 'new grid as' ) ENDIF DO i = 1, layer % idim DO j = 1, layer % jdim IF ( grid % mat (i,j) == grid % nodata ) THEN layer % mat (i,j) = layer % nodata ELSE IF (PRESENT(initial)) THEN layer % mat (i,j) = initial ELSE layer % mat (i,j) = 0. END IF END IF END DO END DO END SUBROUTINE NewGridFloatAsGridFloat !============================================================================== !| Description: ! create a new grid_real using an existing grid_integer as template SUBROUTINE NewGridFloatAsGridInteger & ! (layer, grid, initial) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_integer), INTENT(in) :: grid REAL (KIND = float), OPTIONAL, INTENT(in) :: initial !Arguments with intent(out): TYPE (grid_real), INTENT(OUT) :: layer !!gridreal to return !Local variables: INTEGER (KIND = short) :: ios INTEGER (KIND = short) :: i, j !------------end of declaration------------------------------------------------ layer % jdim = grid % jdim layer % idim = grid % idim layer % varying_mode = 'sequence' !default layer % nodata = MISSING_DEF_REAL layer % valid_min = layer % nodata layer % valid_max = layer % nodata layer % cellsize = grid % cellsize layer % xllcorner = grid % xllcorner layer % yllcorner = grid % yllcorner layer % esri_pe_string = grid % esri_pe_string layer % grid_mapping = grid % grid_mapping layer % reference_time = timeDefault layer % current_time = timeDefault layer % next_time = timeDefault ALLOCATE ( layer % mat ( layer % idim, layer % jdim ), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError, argument = 'new grid as' ) ENDIF DO i = 1, layer % idim DO j = 1, layer % jdim IF ( grid % mat (i,j) == grid % nodata ) THEN layer % mat (i,j) = layer % nodata ELSE IF (PRESENT(initial)) THEN layer % mat (i,j) = initial ELSE layer % mat (i,j) = 0. END IF END IF END DO END DO END SUBROUTINE NewGridFloatAsGridInteger !============================================================================== !| Description: ! create a new grid_integer using an existing grid_integer as template SUBROUTINE NewGridIntegerAsGridInteger & ! (layer, grid, initial) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_integer), INTENT(IN) :: grid INTEGER, OPTIONAL, INTENT(IN) :: initial !Arguments with intent(out): TYPE (grid_integer), INTENT(OUT) :: layer !!grid to be returned !Local variables: INTEGER (KIND = short) :: ios INTEGER (KIND = short) :: i, j !------------end of declaration------------------------------------------------ layer % jdim = grid % jdim layer % idim = grid % idim layer % varying_mode = 'sequence' !default layer % nodata = MISSING_DEF_INT layer % valid_min = layer % nodata layer % valid_max = layer % nodata layer % cellsize = grid % cellsize layer % xllcorner = grid % xllcorner layer % yllcorner = grid % yllcorner layer % esri_pe_string = grid % esri_pe_string layer % grid_mapping = grid % grid_mapping layer % reference_time = timeDefault layer % current_time = timeDefault layer % next_time = timeDefault ALLOCATE ( layer % mat ( layer % idim, layer % jdim ), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError, argument = 'new grid as' ) ENDIF DO i = 1, layer % idim DO j = 1, layer % jdim IF ( grid % mat (i,j) == grid % nodata ) THEN layer % mat (i,j) = layer % nodata ELSE IF (PRESENT(initial)) THEN layer % mat (i,j) = initial ELSE layer % mat (i,j) = 0 END IF END IF END DO END DO END SUBROUTINE NewGridIntegerAsGridInteger !============================================================================== !| Description: ! create a new grid_integer using an existing grid_real as template SUBROUTINE NewGridIntegerAsGridFloat & ! (layer, grid, initial) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_real), INTENT(in) :: grid INTEGER, OPTIONAL, INTENT(in) :: initial !Arguments with intent (out): TYPE (grid_integer), INTENT(OUT) :: layer !!grid to be returned !Local variables: INTEGER (KIND = short) :: ios INTEGER (KIND = short) :: i, j !------------end of declaration------------------------------------------------ layer % jdim = grid % jdim layer % idim = grid % idim layer % varying_mode = 'sequence' !default layer % nodata = MISSING_DEF_INT layer % valid_min = layer % nodata layer % valid_max = layer % nodata layer % cellsize = grid % cellsize layer % xllcorner = grid % xllcorner layer % yllcorner = grid % yllcorner layer % esri_pe_string = grid % esri_pe_string layer % grid_mapping = grid % grid_mapping layer % reference_time = timeDefault layer % current_time = timeDefault layer % next_time = timeDefault ALLOCATE ( layer % mat ( layer % idim, layer % jdim ), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory allocation ', & code = memAllocError, argument = 'new grid as' ) ENDIF DO i = 1, layer % idim DO j = 1, layer % jdim IF ( grid % mat (i,j) == grid % nodata ) THEN layer % mat (i,j) = layer % nodata ELSE IF (PRESENT(initial)) THEN layer % mat (i,j) = initial ELSE layer % mat (i,j) = 0 END IF END IF END DO END DO END SUBROUTINE NewGridIntegerAsGridFloat !============================================================================== !| Description: ! set the standard name of a float grid SUBROUTINE SetStandardNameFloat & ! (name, layer) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: name !Arguments with intent(out): TYPE(grid_real), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % standard_name = name END SUBROUTINE SetStandardNameFloat !============================================================================== !| Description: ! set the standard name of a integer grid SUBROUTINE SetStandardNameInteger & ! (name, layer) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: name !Arguments with intent(out): TYPE(grid_integer), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % standard_name = name END SUBROUTINE SetStandardNameInteger !============================================================================== !| Description: ! set the long name of a float grid SUBROUTINE SetLongNameFloat & ! (name, layer) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: name !Arguments with intent(out): TYPE(grid_real), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % long_name = name END SUBROUTINE SetLongNameFloat !============================================================================== !| Description: ! set the long name of a integer grid SUBROUTINE SetLongNameInteger & ! (name, layer) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: name !Arguments with intent(out): TYPE(grid_integer), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % long_name = name END SUBROUTINE SetLongNameInteger !============================================================================== !| Description: ! set the units of a float grid SUBROUTINE SetUnitsFloat & ! (units, layer) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: units !Arguments with intent(out): TYPE(grid_real), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % units = units END SUBROUTINE SetUnitsFloat !============================================================================== !| Description: ! set the units of a integer grid SUBROUTINE SetUnitsInteger & ! (units, layer) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: units !Arguments with intent(out): TYPE(grid_integer), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % units = units END SUBROUTINE SetUnitsInteger !============================================================================== !| Description: ! set the varying mode of a float grid ! Possible values: ! ! * linear: linear variation between two dates ! * sequence: new grid substitute the previous as frames in a movie ! SUBROUTINE SetVaryingModeFloat & ! (varMod, layer) USE StringManipulation, ONLY: & !Imported routines: StringToLower, StringCompact IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: varMod !Arguments with intent(out): TYPE(grid_real), INTENT(INOUT) :: layer !Local variables CHARACTER (LEN = 20) :: string !------------end of declaration------------------------------------------------ string = StringCompact (StringToLower (varMod) ) IF ( String == 'sequence' .OR. String == 'linear' ) THEN layer % varying_mode = string ELSE CALL Catch ('error', 'GridLib', 'unsupported varying mode: ', & code = unknownOption, argument = string ) END IF END SUBROUTINE SetVaryingModeFloat !============================================================================== !| Description: ! set the varying mode of a integer grid SUBROUTINE SetVaryingModeInteger & ! (varMod, layer) USE StringManipulation, ONLY: & !Imported routines: StringToLower, StringCompact IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: varMod !Arguments with intent(out): TYPE(grid_integer), INTENT(INOUT) :: layer !Local variables CHARACTER (LEN = 20) :: string !------------end of declaration------------------------------------------------ string = StringCompact (StringToLower (varMod) ) IF ( String == 'sequence' .OR. String == 'linear' ) THEN layer % varying_mode = string ELSE CALL Catch ('error', 'GridLib', 'unsupported varying mode: ', & code = unknownOption, argument = string ) END IF END SUBROUTINE SetVaryingModeInteger !============================================================================== !| Description: ! set the reference time of a float grid SUBROUTINE SetReferenceTimeFloat & ! (time, layer) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(IN) :: time !Arguments with intent(out): TYPE(grid_real), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % reference_time = time END SUBROUTINE SetReferenceTimeFloat !============================================================================== !| Description: ! set the reference time of a integer grid SUBROUTINE SetReferenceTimeInteger & ! (time, layer) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(IN) :: time !Arguments with intent(out): TYPE(grid_integer), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % reference_time = time END SUBROUTINE SetReferenceTimeInteger !============================================================================== !| Description: ! set the current time of a float grid SUBROUTINE SetCurrentTimeFloat & ! (time, layer) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(IN) :: time !Arguments with intent(out): TYPE(grid_real), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % current_time = time END SUBROUTINE SetCurrentTimeFloat !============================================================================== !| Description: ! set the current time of a integer grid SUBROUTINE SetCurrentTimeInteger & ! (time, layer) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(IN) :: time !Arguments with intent(out): TYPE(grid_integer), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % current_time = time END SUBROUTINE SetCurrentTimeInteger !============================================================================== !| Description: ! set the esri pe string of a float grid SUBROUTINE SetEsriPeStringFloat & ! (string, layer) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: string !Arguments with intent(out): TYPE(grid_real), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % esri_pe_string = string END SUBROUTINE SetEsriPeStringFloat !============================================================================== !| Description: ! set the esri pe string of a integer grid SUBROUTINE SetEsriPeStringInteger & ! (string, layer) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: string !Arguments with intent(out): TYPE(grid_integer), INTENT(INOUT) :: layer !------------end of declaration------------------------------------------------ layer % esri_pe_string = string END SUBROUTINE SetEsriPeStringInteger !============================================================================== !| Description: ! export grid into netcdf file. Two actions are possible: ! ! * add: add a non-record variable to a netCDF dataset. If file does not exist ! it is created new. The added variable ! must have the same dimensions of already present variables. ! If you try to add a variable already present in netCDF dataset ! an error is returned. ! * append: append more data to record variables along the unlimited dimension. ! If netCDF dataset is created new, dimensions and attributes ! are set, otherwise only new data are written. If current time ! is already present in netcdf dataset an error is returned ! SUBROUTINE ExportGridFloatToNetCDF & ! (grid, file, name, action) USE StringManipulation, ONLY: & ! imported routines: StringToUpper, ToString IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_real), INTENT(IN) :: grid !! grid to be exported CHARACTER (LEN = *), INTENT(IN) :: file !!netcdf file to export to CHARACTER (LEN = *), INTENT(IN) :: name !!name of variable in netcdf CHARACTER (LEN = *), INTENT(IN) :: action !!add or append ! Local variables: INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: ncId !!NetCdf Id for the file INTEGER (KIND = short) :: dimXid !!id of x dimension INTEGER (KIND = short) :: dimYid !!id of y dimension INTEGER (KIND = short) :: dimTid !!id of time dimension INTEGER (KIND = short) :: varXid !!id of x variable INTEGER (KIND = short) :: varYid !!id of y variable INTEGER (KIND = short) :: varTid !!id of time variable INTEGER (KIND = short) :: mapId !!id of grid mapping variable INTEGER (KIND = short), ALLOCATABLE :: varDim (:) !!define dimension of variable INTEGER (KIND = short) :: varId !!id of grid variable CHARACTER (LEN = 25) :: string REAL (KIND = float), ALLOCATABLE :: lats(:), lons(:) !!array containing coordinate INTEGER (KIND = short) :: i !!loop index LOGICAL :: fileExist !!true if netcdf exists LOGICAL :: varExist !!true if variable exists in dataset TYPE (dateTime) :: ref_time !!reference time to calculate time span CHARACTER (LEN = 7) :: dummyString REAL (KIND = float) :: timeSpan !!time between current time and reference time INTEGER (KIND = short) :: records !!number of records stored in netcdf dataset INTEGER (KIND = short) :: slice (3) !!used to write new record in right position INTEGER (KIND = short) :: timeRecord (1) !!to write new time record REAL (KIND = float), ALLOCATABLE :: times (:) !!values stored in time variable LOGICAL :: timeExists REAL (KIND = float), ALLOCATABLE :: tempGrid (:,:) !------------end of declaration------------------------------------------------ !verify if file exists ncStatus = nf90_open (file, NF90_WRITE, ncId) IF (ncStatus == nf90_noerr) THEN fileExist = .TRUE. ELSE fileExist = .FALSE. END IF !if file does not exist create the file IF (.NOT. fileExist) THEN ncStatus = nf90_create ( file, NF90_CLOBBER, ncId ) CALL ncErrorHandler (ncStatus) !define dimensions ncStatus = nf90_def_dim ( ncId, 'x', grid % jdim, dimXid ) CALL ncErrorHandler (ncStatus) ncStatus = nf90_def_dim ( ncId, 'y', grid % idim, dimYid ) CALL ncErrorHandler (ncStatus) !if record (unlimited) variable define time dimension IF ( StringToUpper(action) == 'APPEND' ) THEN ncStatus = nf90_def_dim ( ncId, 'time', NF90_UNLIMITED, dimTid ) CALL ncErrorHandler (ncStatus) END IF !define coordinate variables ncStatus = nf90_def_var (ncId, 'x', NF90_FLOAT, dimXid, varXid) CALL ncErrorHandler (ncStatus) ncStatus = nf90_def_var (ncId, 'y', NF90_FLOAT, dimYid, varYid) CALL ncErrorHandler (ncStatus) !if record (unlimited) variable define time variable IF ( StringToUpper(action) == 'APPEND' ) THEN ncStatus = nf90_def_var ( ncId, 'time', NF90_INT, dimTid, varTid ) CALL ncErrorHandler (ncStatus) END IF !assign attributes to coordinate variables IF (grid % grid_mapping % system == GEODETIC) THEN ncStatus = nf90_put_att (ncId, varXid, 'long_name', 'longitude (centre of grid cell)') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'standard_name', 'longitude') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'units', 'deg') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'axis', 'X') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'long_name', 'latitude (centre of grid cell)') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'standard_name', 'latitude') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'units', 'deg') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'axis', 'Y') CALL ncErrorHandler (ncStatus) ELSE ncStatus = nf90_put_att (ncId, varXid, 'long_name', 'x coordinate of projection') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'standard_name', 'projection_x_coordinate') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'units', 'm') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'axis', 'X') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'long_name', 'y coordinate of projection') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'standard_name', 'projection_y_coordinate') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'units', 'm') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'axis', 'Y') CALL ncErrorHandler (ncStatus) END IF !if record (unlimited) variable define attributes of time variable IF ( StringToUpper(action) == 'APPEND' ) THEN ncStatus = nf90_put_att ( ncId, varTid, 'long_name', 'time' ) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att ( ncId, varTid, 'standard_name', 'time' ) CALL ncErrorHandler (ncStatus) string = grid % current_time ncStatus = nf90_put_att ( ncId, varTid, 'units', 'seconds since ' // string ) CALL ncErrorHandler (ncStatus) END IF !define dummy grid mapping variable ncStatus = nf90_def_var (ncId, 'crs', NF90_INT, mapId) CALL ncErrorHandler (ncStatus) !assign attributes to grid mapping variable ncStatus = nf90_put_att (ncId, mapId, 'grid_mapping_name', & grid % grid_mapping % name) CALL ncErrorHandler (ncStatus) DO i = 1, grid % grid_mapping % basic ncStatus = nf90_put_att (ncId, mapId, grid % grid_mapping % & description (i), grid % grid_mapping % param (i)) CALL ncErrorHandler (ncStatus) END DO !datum ncStatus = nf90_put_att (ncId, mapId, 'datum_code', & grid % grid_mapping % datum % code) CALL ncErrorHandler (ncStatus) !assign global attributes ncStatus = nf90_put_att (ncId, NF90_GLOBAL, 'Conventions', 'CF-1.0') CALL ncErrorHandler (ncStatus) string = UtcNow () ncStatus = nf90_put_att (ncId, NF90_GLOBAL, 'history', 'creation date: ' // string ) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, NF90_GLOBAL, 'institution', 'Politecnico di Milano' ) CALL ncErrorHandler (ncStatus) !end define mode ncStatus = nf90_enddef (ncId) CALL ncErrorHandler (ncStatus) END IF !enter define mode ncStatus = nf90_redef (ncId) CALL ncErrorHandler (ncStatus) !define variable for the grid !If file already exists, have to read dimXd and dimYd !and dimTid if necessary IF ( fileExist) THEN ncStatus = nf90_inq_dimid (ncId, "x", dimXid) CALL ncErrorHandler (ncStatus) ncStatus = nf90_inq_dimid (ncId, "y", dimYid) CALL ncErrorHandler (ncStatus) IF (StringToUpper (action) == 'APPEND') THEN ncStatus = nf90_inq_dimid (ncId, "time", dimTid) CALL ncErrorHandler (ncStatus) ncStatus = nf90_inq_varid (ncId, 'time', varTid) CALL ncErrorHandler (ncStatus) END IF END IF !verify if variable already exists varId = 0 ncStatus = nf90_inq_varid (ncId, name, varId) IF (ncStatus == nf90_noerr) THEN varExist = .TRUE. !if variable already exists in non record dataset: error IF (StringToUpper (action) == 'ADD') THEN CALL Catch ('error', 'GridLib', 'variable ' // TRIM(name) // & ' already exists in netCDf dataset: ', & code = ncIOError, argument = file ) END IF ELSE varExist = .FALSE. !define new variable IF (StringToUpper (action) == 'ADD' ) THEN ALLOCATE ( varDim (2) ) ELSE IF (StringToUpper (action) == 'APPEND' ) THEN ALLOCATE ( varDim (3) ) varDim(3) = dimTid END IF varDim(1) = dimXid varDim(2) = dimYid ncStatus = nf90_def_var (ncId, name, NF90_FLOAT, varDim, varId) CALL ncErrorHandler (ncStatus) DEALLOCATE (varDim) !assign attributes to variable ncStatus = nf90_put_att (ncId, varId, 'long_name', grid % long_name) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varId, 'standard_name', grid % standard_name) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varId, 'units', grid % units) CALL ncErrorHandler (ncStatus) IF (grid % valid_min /= grid % nodata) THEN ncStatus = nf90_put_att (ncId, varId, 'valid_min', grid % valid_min) CALL ncErrorHandler (ncStatus) END IF IF (grid % valid_max /= grid % nodata) THEN ncStatus = nf90_put_att (ncId, varId, 'valid_max', grid % valid_max) CALL ncErrorHandler (ncStatus) END IF ncStatus = nf90_put_att (ncId, varId, '_FillValue', grid % nodata) CALL ncErrorHandler (ncStatus) IF (grid % esri_pe_string /= '') THEN ncStatus = nf90_put_att (ncId, varId, 'esri_pe_string', grid % esri_pe_string) CALL ncErrorHandler (ncStatus) END IF ncStatus = nf90_put_att (ncId, varId, 'varying_mode', grid % varying_mode) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varId, 'grid_mapping', 'crs') CALL ncErrorHandler (ncStatus) END IF !end define mode ncStatus = nf90_enddef (ncId) CALL ncErrorHandler (ncStatus) IF ( .NOT. fileExist) THEN !if file was created new, add coordinate variables !put coordinate ALLOCATE ( lons ( grid % jdim) ) DO i = 1, grid % jdim lons (i) = grid % xllcorner + i * grid % cellsize - grid % cellsize / 2. !lons (i) = int(10000*lons (i) + .5) / 10000.0 END DO ncStatus = nf90_put_var (ncId, varXid, lons ) CALL ncErrorHandler (ncStatus) DEALLOCATE ( lons ) ALLOCATE ( lats ( grid % idim) ) DO i = 1, grid % idim lats (i) = grid % yllcorner + i * grid % cellsize - grid % cellsize / 2. !lats (i) = int(10000*lats (i) + .5) / 10000.0 END DO ncStatus = nf90_put_var (ncId, varYid, lats ) CALL ncErrorHandler (ncStatus) DEALLOCATE ( lats ) END IF !if unlimited variable put time and grid IF ( StringToUpper (action) == 'APPEND' ) THEN !time CALL ParseTime (ncId, ref_time, dummyString) timeSpan = grid % current_time - ref_time !retrieve number of records already stored in dataset ncStatus = nf90_inquire_dimension (ncId, dimTid, len = records) CALL ncErrorHandler (ncStatus) !check if time span already exists ALLOCATE ( times (records) ) ncStatus = nf90_get_var (ncId, varTid, times) timeExists = .FALSE. DO i = 1, records IF ( timeSpan == times (i) ) THEN timeExists = .TRUE. timeRecord = i !(ADD SUPPORT FOR MULTIPLE TIME AXIS !CALL Catch ('error', 'GridLib', & ! 'time is already present in netCDF dataset: ', & ! code = ncIOError, argument = ToString(timeSpan) ) END IF END DO DEALLOCATE (times) !put time span IF ( .NOT. timeExists ) THEN timeRecord = records + 1 ncStatus = nf90_put_var (ncId, varTid, timeSpan, start = timeRecord ) END IF !put grid variable slice = (/ 1, 1, timeRecord /) !swap grid before write to netcdf ALLOCATE (tempGrid (grid % jdim, grid % idim) ) CALL SwapGridRealBack (grid % mat, tempGrid) !ncStatus = nf90_put_var (ncId, varId, tempGrid, start = slice, count = (/1,1,1/) ) ncStatus = nf90_put_var (ncId, varId, tempGrid, start = slice) CALL ncErrorHandler (ncStatus) DEALLOCATE (tempGrid) ELSE IF (StringToUpper (action) == 'ADD' ) THEN !swap grid before write to netcdf ALLOCATE (tempGrid (grid % jdim, grid % idim) ) CALL SwapGridRealBack (grid % mat, tempGrid) !put variable ncStatus = nf90_put_var (ncId, varId, tempGrid ) CALL ncErrorHandler (ncStatus) DEALLOCATE (tempGrid) END IF ! close file ncStatus = nf90_close (ncid) CALL ncErrorHandler (ncStatus) END SUBROUTINE ExportGridFloatToNetCDF !============================================================================== !| Description: ! export grid into netcdf file. Two actions are possible: ! ! * add: add a non-record variable to a netCDF dataset. If file does not exist ! it is created new. The added variable ! must have the same dimensions of already present variables. ! If you try to add a variable already present in netCDF dataset ! an error is returned. ! * append: append more data to record variables along the unlimited dimension. ! If netCDF dataset is created new, dimensions and attributes ! are set, otherwise only new data are written. If current time ! is already present in netcdf dataset an error is returned ! SUBROUTINE ExportGridIntegerToNetCDF & ! (grid, file, name, action) USE StringManipulation, ONLY: & ! imported routines: StringToUpper, ToString IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_integer), INTENT(IN) :: grid !! grid to be exported CHARACTER (LEN = *), INTENT(IN) :: file !!netcdf file to export to CHARACTER (LEN = *), INTENT(IN) :: name !!name of variable in netcdf CHARACTER (LEN = *), INTENT(IN) :: action !!add or append ! Local variables: INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: ncId !!NetCdf Id for the file INTEGER (KIND = short) :: dimXid !!id of x dimension INTEGER (KIND = short) :: dimYid !!id of y dimension INTEGER (KIND = short) :: dimTid !!id of time dimension INTEGER (KIND = short) :: varXid !!id of x variable INTEGER (KIND = short) :: varYid !!id of y variable INTEGER (KIND = short) :: varTid !!id of time variable INTEGER (KIND = short) :: mapId !!id of grid mapping variable INTEGER (KIND = short), ALLOCATABLE :: varDim (:) !!define dimension of variable INTEGER (KIND = short) :: varId !!id of grid variable CHARACTER (LEN = 25) :: string REAL (KIND = float), ALLOCATABLE :: lats(:), lons(:) !!array containing coordinate INTEGER (KIND = short) :: i !!loop index LOGICAL :: fileExist !!true if netcdf exists LOGICAL :: varExist !!true if variable exists in dataset TYPE (dateTime) :: ref_time !!reference time to calculate time span CHARACTER (LEN = 7) :: dummyString REAL (KIND = float) :: timeSpan !!time between current time and reference time INTEGER (KIND = short) :: records !!number of records stored in netcdf dataset INTEGER (KIND = short) :: slice (3) !!used to write new record in right position INTEGER (KIND = short) :: timeRecord (1) !!to write new time record REAL (KIND = float), ALLOCATABLE :: times (:) !!values stored in time variable LOGICAL :: timeExists INTEGER (KIND = long), ALLOCATABLE :: tempGrid (:,:) !------------end of declaration------------------------------------------------ !verify if file exists ncStatus = nf90_open (file, NF90_WRITE, ncId) IF (ncStatus == nf90_noerr) THEN fileExist = .TRUE. ELSE fileExist = .FALSE. END IF !if file does not exist create the file IF (.NOT. fileExist) THEN ncStatus = nf90_create ( file, NF90_CLOBBER, ncId ) CALL ncErrorHandler (ncStatus) !define dimensions ncStatus = nf90_def_dim ( ncId, 'x', grid % jdim, dimXid ) CALL ncErrorHandler (ncStatus) ncStatus = nf90_def_dim ( ncId, 'y', grid % idim, dimYid ) CALL ncErrorHandler (ncStatus) !if record (unlimited) variable define time dimension IF ( StringToUpper(action) == 'APPEND' ) THEN ncStatus = nf90_def_dim ( ncId, 'time', NF90_UNLIMITED, dimTid ) CALL ncErrorHandler (ncStatus) END IF !define coordinate variables ncStatus = nf90_def_var (ncId, 'x', NF90_FLOAT, dimXid, varXid) CALL ncErrorHandler (ncStatus) ncStatus = nf90_def_var (ncId, 'y', NF90_FLOAT, dimYid, varYid) CALL ncErrorHandler (ncStatus) !if record (unlimited) variable define time variable IF ( StringToUpper(action) == 'APPEND' ) THEN ncStatus = nf90_def_var ( ncId, 'time', NF90_FLOAT, dimTid, varTid ) CALL ncErrorHandler (ncStatus) END IF !assign attributes to coordinate variables IF (grid % grid_mapping % system == GEODETIC) THEN ncStatus = nf90_put_att (ncId, varXid, 'long_name', 'longitude (centre of grid cell)') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'standard_name', 'longitude') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'units', 'deg') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'axis', 'X') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'long_name', 'latitude (centre of grid cell)') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'standard_name', 'latitude') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'units', 'deg') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'axis', 'Y') CALL ncErrorHandler (ncStatus) ELSE ncStatus = nf90_put_att (ncId, varXid, 'long_name', 'x coordinate of projection') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'standard_name', 'projection_x_coordinate') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'units', 'm') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varXid, 'axis', 'X') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'long_name', 'y coordinate of projection') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'standard_name', 'projection_y_coordinate') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'units', 'm') CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varYid, 'axis', 'Y') CALL ncErrorHandler (ncStatus) END IF !if record (unlimited) variable define attributes of time variable IF ( StringToUpper(action) == 'APPEND' ) THEN ncStatus = nf90_put_att ( ncId, varTid, 'long_name', 'time' ) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att ( ncId, varTid, 'standard_name', 'time' ) CALL ncErrorHandler (ncStatus) string = grid % current_time ncStatus = nf90_put_att ( ncId, varTid, 'units', 'seconds since ' // string ) CALL ncErrorHandler (ncStatus) END IF !define dummy grid mapping variable ncStatus = nf90_def_var (ncId, 'crs', NF90_INT, mapId) CALL ncErrorHandler (ncStatus) !assign attributes to grid mapping variable ncStatus = nf90_put_att (ncId, mapId, 'grid_mapping_name', & grid % grid_mapping % name) CALL ncErrorHandler (ncStatus) DO i = 1, grid % grid_mapping % basic ncStatus = nf90_put_att (ncId, mapId, grid % grid_mapping % description (i),& grid % grid_mapping % param (i)) CALL ncErrorHandler (ncStatus) END DO !datum ncStatus = nf90_put_att (ncId, mapId, 'datum_code', & grid % grid_mapping % datum % code) CALL ncErrorHandler (ncStatus) !assign global attributes ncStatus = nf90_put_att (ncId, NF90_GLOBAL, 'Conventions', 'CF-1.0') CALL ncErrorHandler (ncStatus) string = UtcNow () ncStatus = nf90_put_att (ncId, NF90_GLOBAL, 'history', 'creation date: ' // string ) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, NF90_GLOBAL, 'institution', 'Politecnico di Milano' ) CALL ncErrorHandler (ncStatus) !end define mode ncStatus = nf90_enddef (ncId) CALL ncErrorHandler (ncStatus) END IF !enter define mode ncStatus = nf90_redef (ncId) CALL ncErrorHandler (ncStatus) !define variable for the grid !If file already exists, have to read dimXd and dimYd !and dimTid if necessary IF ( fileExist) THEN ncStatus = nf90_inq_dimid (ncId, "x", dimXid) CALL ncErrorHandler (ncStatus) ncStatus = nf90_inq_dimid (ncId, "y", dimYid) CALL ncErrorHandler (ncStatus) IF (StringToUpper (action) == 'APPEND') THEN ncStatus = nf90_inq_dimid (ncId, "time", dimTid) CALL ncErrorHandler (ncStatus) ncStatus = nf90_inq_varid (ncId, 'time', varTid) CALL ncErrorHandler (ncStatus) END IF END IF !verify if variable already exists varId = 0 ncStatus = nf90_inq_varid (ncId, name, varId) IF (ncStatus == nf90_noerr) THEN varExist = .TRUE. !if variable already exists in non record dataset: error IF (StringToUpper (action) == 'ADD') THEN CALL Catch ('error', 'GridLib', 'variable ' // TRIM(name) // & ' already exists in netCDf dataset: ', & code = ncIOError, argument = file ) END IF ELSE varExist = .FALSE. !define new variable IF (StringToUpper (action) == 'ADD' ) THEN ALLOCATE ( varDim (2) ) ELSE IF (StringToUpper (action) == 'APPEND' ) THEN ALLOCATE ( varDim (3) ) varDim(3) = dimTid END IF varDim(1) = dimXid varDim(2) = dimYid ncStatus = nf90_def_var (ncId, name, NF90_FLOAT, varDim, varId) CALL ncErrorHandler (ncStatus) DEALLOCATE (varDim) !assign attributes to variable ncStatus = nf90_put_att (ncId, varId, 'long_name', grid % long_name) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varId, 'standard_name', grid % standard_name) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varId, 'units', grid % units) CALL ncErrorHandler (ncStatus) IF (grid % valid_min /= grid % nodata) THEN ncStatus = nf90_put_att (ncId, varId, 'valid_min', grid % valid_min) CALL ncErrorHandler (ncStatus) END IF IF (grid % valid_max /= grid % nodata) THEN ncStatus = nf90_put_att (ncId, varId, 'valid_max', grid % valid_max) CALL ncErrorHandler (ncStatus) END IF ncStatus = nf90_put_att (ncId, varId, '_FillValue', grid % nodata) CALL ncErrorHandler (ncStatus) IF (grid % esri_pe_string /= '') THEN ncStatus = nf90_put_att (ncId, varId, 'esri_pe_string', grid % esri_pe_string) CALL ncErrorHandler (ncStatus) END IF ncStatus = nf90_put_att (ncId, varId, 'varying_mode', grid % varying_mode) CALL ncErrorHandler (ncStatus) ncStatus = nf90_put_att (ncId, varId, 'grid_mapping', 'crs') CALL ncErrorHandler (ncStatus) END IF !end define mode ncStatus = nf90_enddef (ncId) CALL ncErrorHandler (ncStatus) IF ( .NOT. fileExist) THEN !if file was created new, add coordinate variables !put coordinate ALLOCATE ( lons ( grid % jdim) ) DO i = 1, grid % jdim lons (i) = grid % xllcorner + i * grid % cellsize - grid % cellsize / 2. END DO ncStatus = nf90_put_var (ncId, varXid, lons ) CALL ncErrorHandler (ncStatus) DEALLOCATE ( lons ) ALLOCATE ( lats ( grid % idim) ) DO i = 1, grid % idim lats (i) = grid % yllcorner + i * grid % cellsize - grid % cellsize / 2. END DO ncStatus = nf90_put_var (ncId, varYid, lats ) CALL ncErrorHandler (ncStatus) DEALLOCATE ( lats ) END IF !if unlimited variable put time and grid IF ( StringToUpper (action) == 'APPEND' ) THEN !time CALL ParseTime (ncId, ref_time, dummyString) timeSpan = grid % current_time - ref_time !retrieve number of records already stored in dataset ncStatus = nf90_inquire_dimension (ncId, dimTid, len = records) CALL ncErrorHandler (ncStatus) !check if time span already exists ALLOCATE ( times (records) ) ncStatus = nf90_get_var (ncId, varTid, times) timeExists = .FALSE. DO i = 1, records IF ( timeSpan == times (i) ) THEN timeExists = .TRUE. timeRecord = i !(ADD SUPPORT FOR MULTIPLE TIME AXIS !CALL Catch ('error', 'GridLib', & ! 'time is already present in netCDF dataset: ', & ! code = ncIOError, argument = ToString(timeSpan) ) END IF END DO DEALLOCATE (times) !put time span IF ( .NOT. timeExists ) THEN timeRecord = records + 1 ncStatus = nf90_put_var (ncId, varTid, timeSpan, start = timeRecord ) END IF !put grid variable slice = (/ 1, 1, timeRecord /) !swap grid before write to netcdf ALLOCATE (tempGrid (grid % jdim, grid % idim) ) CALL SwapGridIntegerBack (grid % mat, tempGrid) ncStatus = nf90_put_var (ncId, varId, tempGrid, start = slice ) CALL ncErrorHandler (ncStatus) DEALLOCATE (tempGrid) ELSE IF (StringToUpper (action) == 'ADD' ) THEN !swap grid before write to netcdf ALLOCATE (tempGrid (grid % jdim, grid % idim) ) CALL SwapGridIntegerBack (grid % mat, tempGrid) !put variable ncStatus = nf90_put_var (ncId, varId, tempGrid ) CALL ncErrorHandler (ncStatus) DEALLOCATE (tempGrid) END IF ! close file ncStatus = nf90_close (ncid) CALL ncErrorHandler (ncStatus) END SUBROUTINE ExportGridIntegerToNetCDF !============================================================================== !| Description: ! export grid_real to file. ! List of supported format: ! ! * `ESRI_ASCII`: ESRI ASCII GRID ! * `ESRI_BINARY`: ESRI BINARY GRID ! SUBROUTINE ExportGridFloatToFile & ! (layer, fileName, fileFormat) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_real), INTENT (IN) :: layer CHARACTER (LEN = *), INTENT (IN) :: fileName INTEGER (KIND = short), INTENT(IN) :: fileFormat !------------end of declaration------------------------------------------------ IF ( fileformat == ESRI_ASCII ) THEN CALL ExportGridFloatToESRI_ASCII (layer, fileName) ELSE IF ( fileformat == ESRI_BINARY ) THEN CALL ExportGridFloatToESRI_BINARY (layer, fileName) ELSE CALL Catch ('error', 'GridLib', & 'unknown option on exporting file: ', & code = unknownOption, argument = fileName ) END IF END SUBROUTINE ExportGridFloatToFile !============================================================================== !| Description: ! export grid_real to a ESRI ASCII file. SUBROUTINE ExportGridFloatToESRI_ASCII & ! (layer, fileName) USE Utilities, ONLY: & !Imported routines: GetUnit IMPLICIT NONE !Arguments with intent(in): TYPE (grid_real), INTENT (IN) :: layer CHARACTER (LEN = *), INTENT (IN) :: fileName !Local variables: INTEGER (KIND = short) :: fileUnit INTEGER (KIND = short) :: ios INTEGER (KIND = short) :: i,j !------------end of declaration------------------------------------------------ !open file fileUnit = GetUnit () OPEN (UNIT = fileUnit, file = fileName, IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = fileName ) END IF !write header WRITE(fileUnit,'(a14,i10)') "ncols ", layer % jdim WRITE(fileUnit,'(a14,i10)') "nrows ", layer % idim WRITE(fileUnit,'(a14,f15.5)') "xllcorner ", layer % xllcorner WRITE(fileUnit,'(a14,f15.5)') "yllcorner ", layer % yllcorner WRITE(fileUnit,'(a14,f15.5)') "cellsize ", layer % cellsize WRITE(fileUnit,'(a14,E14.7)') "NODATA_value ", layer % nodata !write data DO i = 1, layer % idim DO j = 1, layer % jdim - 1 WRITE(fileUnit,'(E14.7," ")', ADVANCE = 'no') layer % mat(i,j) END DO WRITE(fileUnit,'(E14.7," ")') layer % mat(i,layer % jdim) END DO CLOSE (fileUnit) END SUBROUTINE ExportGridFloatToESRI_ASCII !============================================================================== !| Description: ! export grid_real to a ESRI BINARY file. SUBROUTINE ExportGridFloatToESRI_BINARY & ! (layer, fileName) USE Utilities, ONLY: & !Imported routines: GetUnit IMPLICIT NONE !Arguments with intent(in): TYPE (grid_real), INTENT (IN) :: layer CHARACTER (LEN = *), INTENT (IN) :: fileName !Local variables: INTEGER (KIND = short) :: fileUnit INTEGER (KIND = short) :: ios INTEGER (KIND = short) :: recordLength INTEGER (KIND = short) :: recordNumber INTEGER (KIND = short) :: i,j !------------end of declaration------------------------------------------------ !open file fileUnit = GetUnit () OPEN (UNIT = fileUnit, file = TRIM(fileName) // '.hdr', IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = TRIM(fileName) // '.hdr' ) END IF !write header WRITE(fileUnit,'(a14,i10)') "ncols ", layer % jdim WRITE(fileUnit,'(a14,i10)') "nrows ", layer % idim WRITE(fileUnit,'(a14,f15.5)') "xllcorner ", layer % xllcorner WRITE(fileUnit,'(a14,f15.5)') "yllcorner ", layer % yllcorner WRITE(fileUnit,'(a14,f15.5)') "cellsize ", layer % cellsize WRITE(fileUnit,'(a14,E14.7)') "NODATA_value ", layer % nodata WRITE(fileUnit,'(a18)') "BYTEORDER LSBFIRST" CLOSE (fileUnit) !write data fileUnit = GetUnit () INQUIRE (IOLENGTH = recordLength) 100_float OPEN (UNIT = fileUnit, file = TRIM(fileName) // '.flt', & FORM = 'UNFORMATTED', ACCESS = 'DIRECT', RECL = recordLength, & IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = TRIM(fileName) // '.flt' ) END IF recordNumber = 0 DO i = 1,layer % idim DO j = 1, layer % jdim recordNumber = recordNumber + 1 WRITE (fileUnit,REC = recordNumber) layer % mat (i,j) END DO END DO CLOSE (fileUnit) END SUBROUTINE ExportGridFloatToESRI_BINARY !============================================================================== !| Description: ! export grid_integer to file. ! List of supported format: ! ! * `ESRI_ASCII`: ESRI ASCII GRID ! * `ESRI_BINARY`: ESRI BINARY GRID ! SUBROUTINE ExportGridIntegerToFile & ! (layer, fileName, fileFormat) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_integer), INTENT (IN) :: layer CHARACTER (LEN = *), INTENT (IN) :: fileName INTEGER (KIND = short), INTENT(IN) :: fileFormat !------------end of declaration------------------------------------------------ IF ( fileformat == ESRI_ASCII ) THEN CALL ExportGridIntegerToESRI_ASCII (layer, fileName) ELSE IF ( fileformat == ESRI_BINARY ) THEN CALL ExportGridIntegerToESRI_BINARY (layer, fileName) ELSE CALL Catch ('error', 'GridLib', & 'unknown option on exporting file: ', & code = unknownOption, argument = fileName ) END IF END SUBROUTINE ExportGridIntegerToFile !============================================================================== !| Description: ! export grid_integer to a ESRI ASCII file. SUBROUTINE ExportGridIntegerToESRI_ASCII & ! (layer, fileName) USE Utilities, ONLY: & !Imported routines: GetUnit IMPLICIT NONE !Arguments with intent(in): TYPE (grid_integer), INTENT (IN) :: layer CHARACTER (LEN = *), INTENT (IN) :: fileName !Local variables: INTEGER (KIND = short) :: fileUnit INTEGER (KIND = short) :: ios INTEGER (KIND = short) :: i,j !------------end of declaration------------------------------------------------ !open file fileUnit = GetUnit () OPEN (UNIT = fileUnit, file = fileName, IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = fileName ) END IF !write header WRITE(fileUnit,'(a14,i4)') "ncols ", layer % jdim WRITE(fileUnit,'(a14,i4)') "nrows ", layer % idim WRITE(fileUnit,'(a14,f15.5)') "xllcorner ", layer % xllcorner WRITE(fileUnit,'(a14,f15.5)') "yllcorner ", layer % yllcorner WRITE(fileUnit,'(a14,f15.5)') "cellsize ", layer % cellsize WRITE(fileUnit,'(a14,i8)') "NODATA_value ", layer % nodata !write data DO i = 1,layer % idim DO j = 1, layer % jdim - 1 WRITE(fileUnit,'(i8," ")', ADVANCE = 'no') layer % mat(i,j) END DO WRITE(fileUnit,'(i8," ")') layer % mat(i,layer % jdim) END DO CLOSE (fileUnit) END SUBROUTINE ExportGridIntegerToESRI_ASCII !============================================================================== !| Description: ! export grid_integer to a ESRI BINARY file. SUBROUTINE ExportGridIntegerToESRI_BINARY & ! (layer, fileName) USE Utilities, ONLY: & !Imported routines: GetUnit IMPLICIT NONE !Arguments with intent(in): TYPE (grid_integer), INTENT (IN) :: layer CHARACTER (LEN = *), INTENT (IN) :: fileName !Local variables: INTEGER (KIND = short) :: fileUnit INTEGER (KIND = short) :: ios INTEGER (KIND = short) :: recordLength INTEGER (KIND = short) :: recordNumber INTEGER (KIND = short) :: i,j !------------end of declaration------------------------------------------------ !open file fileUnit = GetUnit () OPEN (UNIT = fileUnit, file = TRIM(fileName) // '.hdr', IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = TRIM(fileName) // '.hdr' ) END IF !write header WRITE(fileUnit,'(a14,i4)') "ncols ", layer % jdim WRITE(fileUnit,'(a14,i4)') "nrows ", layer % idim WRITE(fileUnit,'(a14,f15.5)') "xllcorner ", layer % xllcorner WRITE(fileUnit,'(a14,f15.5)') "yllcorner ", layer % yllcorner WRITE(fileUnit,'(a14,f15.5)') "cellsize ", layer % cellsize WRITE(fileUnit,'(a14,i8)') "NODATA_value ", layer % nodata WRITE(fileUnit,'(a18)') "BYTEORDER LSBFIRST" CLOSE (fileUnit) !write data fileUnit = GetUnit () INQUIRE (IOLENGTH = recordLength) 100_long OPEN (UNIT = fileUnit, file = TRIM(fileName) // '.flt', & FORM = 'UNFORMATTED', ACCESS = 'DIRECT', RECL = recordLength, & IOSTAT = ios) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'error opening file: ', & code = openFileError, argument = TRIM(fileName) // '.flt' ) END IF recordNumber = 0 DO i = 1,layer % idim DO j = 1, layer % jdim recordNumber = recordNumber + 1 WRITE (fileUnit,REC = recordNumber) layer % mat (i,j) END DO END DO CLOSE (fileUnit) END SUBROUTINE ExportGridIntegerToESRI_BINARY !============================================================================== !| Description: ! read and calculate georeferencing informations from netCDF file. SUBROUTINE GetGeoreferenceFromNCdataSet & ! (ncId, varId, cellsize, xll, yll, grid_mapping, Iraster, Jraster, point2raster) USE StringManipulation, ONLY: & !Imported routines: StringCompact, & StringToLower IMPLICIT NONE ! Arguments with intent (in): INTEGER (KIND = short), INTENT(IN) :: ncId !!NetCdf Id for the file INTEGER (KIND = short), INTENT(IN) :: varId !!variable Id LOGICAL, INTENT(IN) :: point2raster ! Arguments with intent(out): REAL (KIND = float), INTENT(OUT) :: cellsize !!cell resolution REAL (KIND = float), INTENT(OUT) :: xll !!x coordinate of lower left corner of map REAL (KIND = float), INTENT(OUT) :: yll !!y coordinate of lower left corner of map TYPE (CRS), INTENT(OUT) :: grid_mapping !!coordinate reference system ! Arguments with intent(inout): INTEGER (KIND = short), INTENT(INOUT) :: Iraster (:,:) INTEGER (KIND = short), INTENT(INOUT) :: Jraster (:,:) ! Local variables: INTEGER (KIND = short) :: x, y !!number of columns and rows INTEGER (KIND = short) :: idx, idy REAL (KIND = float), ALLOCATABLE :: xArray (:), yArray (:) REAL (KIND = float), ALLOCATABLE :: xMatrix (:,:), yMatrix (:,:) REAL (KIND = float), ALLOCATABLE :: xMatrix2 (:,:), yMatrix2 (:,:) REAL (KIND = float) :: jdim, idim !! x and y resolution REAL (KIND = float) :: xmin !! minimum value of x coordinate in xMatrix REAL (KIND = float) :: xmax !! maximum value of x coordinate in xMatrix REAL (KIND = float) :: ymin !! minimum value of y coordinate in yMatrix REAL (KIND = float) :: ymax !! maximum value of y coordinate in yMatrix INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: i, j INTEGER (KIND = short) :: r, c CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 100) :: variableName CHARACTER (LEN = 1) :: shp !!define if x and y are vectors or matrix INTEGER (KIND = short) :: mappingId !!grid mapping variable Id INTEGER (KIND = short) :: refSystem !!reference system INTEGER (KIND = short) :: datum !!datum code REAL (KIND = float) :: lat0 !!latitude of projection origin REAL (KIND = float) :: centM !!longitude of central meridian REAL (KIND = float) :: primeM !!longitude of prime meridian REAL (KIND = float) :: falseE REAL (KIND = float) :: falseN REAL (KIND = float) :: scale_factor REAL (KIND = float) :: tol = 0.01 !!tolerance for checking regularly spaced grid INTEGER (KIND = short) :: latlon !------------end of declaration------------------------------------------------ !get length of x and y sizes and related id CALL GetXYSizesFromFile (ncId, x, y, idx = idx, idy = idy, shape = shp, latlon = latlon) IF (shp == 'v' ) THEN !x and y values are stored in vectors !allocate space to read coordinates ALLOCATE ( xArray (x) ) ALLOCATE ( yArray (y) ) !read x coordinates ncStatus = nf90_get_var (ncId, idx, xArray) CALL ncErrorHandler (ncStatus) !read y coordinates ncStatus = nf90_get_var (ncId, idy, yArray) CALL ncErrorHandler (ncStatus) !calculate and check x spatial resolution jdim = xArray (2) - xArray (1) DO i = 3, x IF ( ABS((xArray (i) - xArray (i-1) - jdim ) / jdim) > tol) THEN CALL Catch ('error', 'GridLib', 'x dimension not regularly spaced' ) END IF END DO !calculate and check y spatial resolution idim = yArray (2) - yArray (1) DO i = 3, y IF ( ABS((yArray (i) - yArray (i-1) - idim) / idim) > tol) THEN CALL Catch ('error', 'GridLib', 'y dimension not regularly spaced' ) END IF END DO !set cellsize IF (ABS((jdim - idim)/jdim) < tol) THEN cellsize = (jdim + idim) / 2. ELSE CALL Catch ('error', 'GridLib', 'x resolution different from y resolution' ) END IF !calculate xll and yll xll = MINVAL (xArray) - cellsize / 2. yll = MINVAL (yArray) - cellsize / 2. !deallocate arrays DEALLOCATE ( xArray ) DEALLOCATE ( yArray ) ELSE IF (shp == 'm') THEN !x and y values are stored in matrix !allocate space to read coordinates ALLOCATE ( xMatrix (x,y) ) ALLOCATE ( yMatrix (x,y) ) ALLOCATE ( xMatrix2 (y,x) ) ALLOCATE ( yMatrix2 (y,x) ) !read x coordinates ncStatus = nf90_get_var (ncId, idx, xMatrix) CALL ncErrorHandler (ncStatus) !read y coordinates ncStatus = nf90_get_var (ncId, idy, yMatrix) CALL ncErrorHandler (ncStatus) !calculate x spatial resolution xmin = MINVAL(xMatrix) xmax = MAXVAL(xMatrix) jdim = (xmax - xmin) / (x - 1) !calculate y spatial resolution ymin = MINVAL(yMatrix) ymax = MAXVAL(yMatrix) idim = (ymax - ymin) / (y - 1) !set cellsize !cellsize = (idim + jdim) / 2. cellsize = MAX(idim , jdim) !calculate xll and yll xll = xmin - cellsize / 2. yll = ymin - cellsize / 2. !swap x and y matrix CALL SwapGridRealForward (xMatrix, xMatrix2, latlon) CALL SwapGridRealForward (yMatrix, yMatrix2, latlon) !build Iraster and Jraster matrix IF (point2raster) THEN Iraster = -9999 Jraster = -9999 DO i = 1, y DO j = 1, x c = INT ( (xMatrix2(i,j) - xll) / cellsize ) + 1 r = INT ( (yll + y * cellsize - yMatrix2(i,j)) / cellsize ) + 1 IF (r >= 1 .AND. r <= y .AND. c >= 1 .AND. c <= x ) THEN Iraster(i,j) = r Jraster(i,j) = c END IF END DO END DO END IF !deallocate DEALLOCATE (xMatrix) DEALLOCATE (yMatrix) DEALLOCATE (xMatrix2) DEALLOCATE (yMatrix2) END IF !read coordinate reference system attribute = '' ncStatus = nf90_get_att (ncId, varid = varId, name = 'grid_mapping', & values = attribute) IF (ncStatus == nf90_noerr) THEN !check if grid_mapping exists !retrieve varId of grid_mapping variable ncStatus = nf90_inq_varid (ncId, attribute, mappingId) CALL ncErrorHandler (ncStatus) !retrieve datum EPSG code (extension to CF conventions) ncStatus = nf90_get_att (ncId, varid = mappingId, name = 'datum_code', & values = datum) !if datum is not specified, datum_code does not exist, set default to WGS84 IF (ncStatus /= nf90_noerr) THEN datum = WGS84 END IF !retrieve grid mapping name ncStatus = nf90_get_att (ncId, varid = mappingId, name = 'grid_mapping_name', & values = attribute) IF (ncStatus /= nf90_noerr) THEN !grid_mapping_name not found !set default to geodetic WGS84 CALL SetCRS (GEODETIC, WGS84, grid_mapping) CALL SetGeodeticParameters (grid_mapping, prime_meridian = 0.0) ELSE !retrieve reference system parameters IF ( StringToLower (StringCompact (attribute) ) == 'transverse_mercator' ) THEN CALL SetCRS (TM, datum, grid_mapping ) ncStatus = nf90_get_att (ncId, varid = mappingId, & name = 'longitude_of_central_meridian', & values = centM) ncStatus = nf90_get_att (ncId, varid = mappingId, & name = 'central_meridian', & values = centM) ncStatus = nf90_get_att (ncId, varid = mappingId, & name = 'latitude_of_projection_origin', & values = lat0) ncStatus = nf90_get_att (ncId, varid = mappingId, & name = 'scale_factor_at_central_meridian', & values = scale_factor) ncStatus = nf90_get_att (ncId, varid = mappingId, & name = 'scale_factor', & values = scale_factor) ncStatus = nf90_get_att (ncId, varid = mappingId, & name = 'false_easting', & values = falseE) ncStatus = nf90_get_att (ncId, varid = mappingId, & name = 'false_northing', & values = falseN) CALL SetTransverseMercatorParameters (grid_mapping, lat0, centM, & falseE, falseN, scale_factor) ELSE IF ( StringCompact (attribute) == 'latitude_longitude' ) THEN CALL SetCRS (GEODETIC, datum, grid_mapping ) ncStatus = nf90_get_att (ncId, varid = mappingId, & name = 'longitude_of_prime_meridian', & values = primeM) CALL SetGeodeticParameters (grid_mapping, primeM) ELSE !other reference systems not yet supported END IF END IF ELSE !CRS not specified: set default to geodetic WGS84 !retrieve name of variable ncStatus = nf90_inquire_variable(ncId, varId, name = variableName) CALL ncErrorHandler (ncStatus) CALL Catch ('warning', 'GridLib', & 'grid_mapping not found in variable: ', & argument = variableName ) CALL Catch ('info', 'GridLib', & 'set default coordinate reference system to geodetic wgs84') CALL SetCRS (GEODETIC, WGS84, grid_mapping) CALL SetGeodeticParameters (grid_mapping, prime_meridian = 0.0) END IF RETURN END SUBROUTINE GetGeoreferenceFromNCdataSet !============================================================================== !| Description: ! extraxt the number of columns (x) and rows (y) from netCDF file. SUBROUTINE GetXYSizesFromFile & ! (ncId, x, y, idx, idy, shape, latlon) IMPLICIT NONE !Arguments with intent(in): INTEGER (KIND = short), INTENT(IN) :: ncId !!NetCdf Id for the file !Arguments with intent(out): INTEGER (KIND = short), INTENT(OUT) :: x, y !!number of columns and rows INTEGER (KIND = short), OPTIONAL, INTENT(OUT) :: idx !! id of x variable INTEGER (KIND = short), OPTIONAL, INTENT(OUT) :: idy !! id of y variable CHARACTER (LEN = *), INTENT(OUT) :: shape !!ccordinate in vector or matrix INTEGER (KIND = short), OPTIONAL, INTENT(OUT) :: latlon !! coordinate order is lat-lon (1) or lon-lat (2) ! Local variables: INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: IdXCoordinate !!Id of the variable containing !!information on x ccordinate INTEGER (KIND = short) :: IdYCoordinate !!Id of the variable containing !!information on y ccordinate INTEGER (KIND = short) :: IdXDim !!Id of the X dimension INTEGER (KIND = short) :: IdYDim !!Id of the Y dimension INTEGER (KIND = short) :: nDimsVar !!number of dimensions of a variable INTEGER (KIND = short), ALLOCATABLE :: dimIds (:) !!dimension Ids of a variable CHARACTER (LEN = 80) :: attribute INTEGER (KIND = short) :: i LOGICAL :: stdnamefound LOGICAL :: longnamefound LOGICAL :: unitsfound !------------end of declaration------------------------------------------------ !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nVariables = nVars ) CALL ncErrorHandler (ncStatus) stdnamefound = .FALSE. longnamefound = .FALSE. unitsfound = .FALSE. !scan dataset searching for variable referring to grid mapping DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN !'standard_name' is found IF ( attribute(1:23) == 'projection_x_coordinate' .OR. & !projected reference system attribute(1:9) == 'longitude' .OR. & !geographic reference system attribute(1:7) == 'Easting' .OR. & attribute(1:14) == 'grid_longitude')THEN !rotated pole system stdnamefound = .TRUE. IdXCoordinate = i ELSE IF ( attribute(1:23) == 'projection_y_coordinate' .OR. & attribute(1:8) == 'latitude' .OR. & !geographic reference system attribute(1:8) == 'Northing' .OR. & attribute(1:13) == 'grid_latitude')THEN !rotated pole system stdnamefound = .TRUE. IdYCoordinate = i END IF END IF END DO IF ( .NOT. stdnamefound) THEN CALL Catch ('warning', 'GridLib', & 'standard name not found while searching for x and y axis') END IF DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'long_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN !'long_name' is found IF ( attribute(1:23) == 'projection_x_coordinate' .OR. & !projected reference system attribute(1:9) == 'longitude' .OR. & !geographic reference system attribute(1:7) == 'Easting' .OR. & attribute(1:14) == 'grid_longitude')THEN !rotated pole system longnamefound = .TRUE. IdXCoordinate = i ELSE IF ( attribute(1:23) == 'projection_y_coordinate' .OR. & attribute(1:8) == 'latitude' .OR. & !geographic reference system attribute(1:8) == 'Northing' .OR. & attribute(1:13) == 'grid_latitude')THEN !rotated pole system longnamefound = .TRUE. IdYCoordinate = i END IF END IF END DO IF ( .NOT. longnamefound) THEN CALL Catch ('warning', 'GridLib', & 'long_name not found while searching for x and y axis') END IF DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'units', & values = attribute) IF (ncStatus == nf90_noerr) THEN !'units' is found IF ( attribute(1:11) == 'degree_east' )THEN unitsfound = .TRUE. IdXCoordinate = i ELSE IF ( attribute(1:12) == 'degree_north' )THEN unitsfound = .TRUE. IdYCoordinate = i END IF END IF END DO IF ( .NOT. unitsfound) THEN CALL Catch ('warning', 'GridLib', & 'units not found while searching for x and y axis') END IF IF ( .NOT. stdnamefound .AND. & .NOT. longnamefound .AND. & .NOT. unitsfound ) THEN CALL Catch ('error', 'GridLib', & 'x and y axis can not be found') END IF IF (PRESENT (idx) ) THEN idx = IdXCoordinate END IF IF (PRESENT (idy) ) THEN idy = IdYCoordinate END IF !retrieve Id of X dimension ncStatus = nf90_inquire_variable (ncId, IdXCoordinate, ndims = nDimsVar) CALL ncErrorHandler (ncStatus) ALLOCATE (dimIds (nDimsVar)) ncStatus = nf90_inquire_variable (ncId, IdXCoordinate, dimids = dimIds) CALL ncErrorHandler (ncStatus) IF ( nDimsVar == 1) THEN !x is a vector IdXDim = dimIds (1) shape = 'v' ELSE IF ( nDimsVar == 2) THEN !x is a matrix. Used for not regular grids. !Suppose that x is the first dimension IdXDim = dimIds (1) shape = 'm' END IF DEALLOCATE (dimIds) !retrieve Id of Y dimension ncStatus = nf90_inquire_variable (ncId, IdYCoordinate, ndims = nDimsVar) CALL ncErrorHandler (ncStatus) ALLOCATE (dimIds (nDimsVar)) ncStatus = nf90_inquire_variable (ncId, IdYCoordinate, dimids = dimIds) CALL ncErrorHandler (ncStatus) IF ( nDimsVar == 1) THEN !y is a vector IdYDim = dimIds (1) ELSE IF ( nDimsVar == 2) THEN !y is a matrix. Used for not regular grids. !Suppose that y is the second dimension IdYDim = dimIds (2) END IF DEALLOCATE (dimIds) !set coordinate order convention lon-lat or lat-lon IF (PRESENT (latlon) ) THEN IF (IdXDim > IdYDim) THEN !coordinate convention order is lon-lat latlon = 1 !latlon = 2 !DAO ELSE !coordinate convention order is lat-lon latlon = 2 !latlon = 1 !DAO END IF END IF !retrieve length ncStatus = nf90_inquire_dimension (ncId, IdXDim, len = x) CALL ncErrorHandler (ncStatus) ncStatus = nf90_inquire_dimension (ncId, IdYDim, len = y) CALL ncErrorHandler (ncStatus) RETURN END SUBROUTINE GetXYSizesFromFile !============================================================================== !| Description: ! Calculate the index to extract the corresponding slice from netcdf file. FUNCTION TimeIndex & ! (ncId, refTime, timeUnit, time) & RESULT (index) USE Units, ONLY: & ! Imported parameters: minute, hour, day, month USE StringManipulation, ONLY: & !Imported routines: StringtOsHORT IMPLICIT NONE ! Arguments with intent(in): INTEGER (KIND = short), INTENT(IN) :: ncId !!NetCdf Id for the file TYPE (DateTime), INTENT(IN) :: refTime !!reference time to calculate time index CHARACTER (LEN = *) :: timeUnit TYPE (DateTime), INTENT (IN) :: time !!time to calculate index from reference time ! Local variables: INTEGER (KIND = short) :: index INTEGER (KIND = long) :: difference INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: nDims !!number of dimensions INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: nAtts !!number of global attributes INTEGER (KIND = short) :: length !!length of time dimension INTEGER (KIND = short) :: idTime !!Id of the variable containing !!information on time ccordinate INTEGER (KIND = short) :: timeNumber INTEGER (KIND = short) :: CurrentTimeNumber CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 100) :: variableName INTEGER (KIND = short) :: i !!loop index INTEGER :: slice (2) INTEGER :: current LOGICAL :: found CHARACTER (LEN = 25) :: string CHARACTER (LEN = 19) :: str INTEGER, DIMENSION(NF90_MAX_VAR_DIMS) :: dimIDs !------------end of declaration------------------------------------------------ !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nDimensions = nDims, & nVariables = nVars, & nAttributes = nAtts ) CALL ncErrorHandler (ncStatus) !search for time variable DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN IF ( attribute(1:4) == 'time' ) THEN idTime = i EXIT END IF ELSE !standard_name is not defined: search for variable named 'time' !ncStatus = nf90_inq_varid (ncId, 'time', varid = i ) ncstatus = nf90_inquire_variable(ncId, varId = i, name = variableName) IF (LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'time' .OR. & LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'Time' .OR. & LEN_TRIM(variableName) == 5 .AND. & variableName(1:5) == 'Times' ) THEN !variable 'time' found idTime = i EXIT END IF END IF END DO !retrieve time length length = GetTimeSteps (ncId) !ncStatus = nf90_inquire_variable(ncid, idTime, dimids = dimIDs) !CALL ncErrorHandler (ncStatus) ! !ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(2), len = length) !CALL ncErrorHandler (ncStatus) !search for current time found = .FALSE. IF (DateTimeIsDefault(refTime)) THEN !build datetime string as used in netcdf file i.e 2007-10-11_00:00:00 timeString = time !2000-01-01T00:00:00+00:00 timeString = timeString(1:10) // '_' // timeString(12:19) !CurrentTimeNumber = StringToShort (timeString) DO i = 1, length slice(1) = 1 slice(2) = i ncStatus = nf90_get_var (ncId, idTime, str , start = slice) CALL ncErrorHandler (ncStatus) IF (TRIM(str) == TRIM(timeString)) THEN found = .TRUE. index = i EXIT END IF END DO ELSE !calculate time span in appropriate unit difference = time - refTime SELECT CASE (timeUnit) CASE ('minutes') difference = difference / INT(minute) CASE ('hours') difference = difference / INT(hour) CASE ('days') difference = difference / INT(day) CASE ('months') difference = difference / INT(month) END SELECT DO i = 1, length slice(1) = i ncStatus = nf90_get_var (ncId, idTime, current , start = slice) CALL ncErrorHandler (ncStatus) IF (current == difference) THEN found = .TRUE. index = i EXIT END IF END DO END IF IF ( .NOT. found ) THEN string = time CALL Catch ('error', 'GridLib', & 'time not found in netcdf file: ', & argument = string ) END IF END FUNCTION TimeIndex !============================================================================== !| Description ! returns the time of the next grid in netcdf dataset FUNCTION NextTime & ! (ncId, refTime, timeUnit, current) & ! RESULT (time) USE Units, ONLY: & ! Imported parameters: minute, hour, day, month USE StringManipulation, ONLY: & !Imported routines: ToString IMPLICIT NONE ! Arguments with intent(in): INTEGER (KIND = short), INTENT(IN) :: ncId !!NetCdf Id for the file TYPE (DateTime), INTENT(IN) :: refTime !!reference time to calculate time index CHARACTER (LEN = *), INTENT(IN) :: timeUnit INTEGER (KIND = short), INTENT(IN) :: current !!current time step !Local variables: TYPE (DateTime) :: time !!returned time of the next grid INTEGER (KIND = short) :: next !!index of next time step INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: nDims !!number of dimensions INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: nAtts !!number of global attributes INTEGER (KIND = short) :: length !!length of time dimension INTEGER (KIND = short) :: idTime !!Id of the variable containing !!information on time ccordinate CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 100) :: variableName INTEGER (KIND = short) :: i !!loop index INTEGER :: slice (2) INTEGER :: timeSpan CHARACTER (LEN = 80) :: string CHARACTER (LEN = 19) :: str INTEGER, DIMENSION(NF90_MAX_VAR_DIMS) :: dimIDs !------------end of declaration------------------------------------------------ !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nDimensions = nDims, & nVariables = nVars, & nAttributes = nAtts ) CALL ncErrorHandler (ncStatus) !search for time variable DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN IF ( attribute(1:4) == 'time' ) THEN idTime = i EXIT END IF ELSE !standard_name is not defined: search for variable named 'time' !ncStatus = nf90_inq_varid (ncId, 'time', varid = i ) ncstatus = nf90_inquire_variable(ncId, varId = i, name = variableName) IF (LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'time' .OR. & LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'Time' .OR. & LEN_TRIM(variableName) == 5 .AND. & variableName(1:5) == 'Times') THEN !variable 'time' found idTime = i EXIT END IF END IF END DO !inquire time length ncStatus = nf90_inquire_variable(ncid, idTime, dimids = dimIDs) CALL ncErrorHandler (ncStatus) !ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(2), len = length) ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(1), len = length) CALL ncErrorHandler (ncStatus) !set next time step IF (current < length) THEN next = current + 1 ELSE next = length END IF !compute date corresponding to next time step !slice(1) = 1 !slice(2) = next slice(1) = next slice(2) = 1 IF (DateTimeIsDefault(refTime)) THEN ncStatus = nf90_get_var (ncId, idTime, str , start = slice) CALL ncErrorHandler (ncStatus) string = str(1:10) // 'T' // str(12:19) // '+00:00' time = string ELSE ncStatus = nf90_get_var (ncId, idTime, timeSpan , start = slice) CALL ncErrorHandler (ncStatus) SELECT CASE (timeUnit) CASE ('minutes') timeSpan = timeSpan * minute CASE ('hours') timeSpan = timeSpan * hour CASE ('days') timeSpan = timeSpan * day CASE ('months') timeSpan = timeSpan * month END SELECT time = refTime + timeSpan END IF RETURN END FUNCTION NextTime !============================================================================== !| Description ! returns time given time index FUNCTION TimeByIndex & ! (ncId, refTime, timeUnit, index) & ! RESULT (time) USE Units, ONLY: & ! Imported parameters: minute, hour, day, month USE StringManipulation, ONLY: & !Imported routines: ToString IMPLICIT NONE ! Arguments with intent(in): INTEGER (KIND = short), INTENT(IN) :: ncId !!NetCdf Id for the file TYPE (DateTime), INTENT(IN) :: refTime !!reference time to calculate time index CHARACTER (LEN = *), INTENT(IN) :: timeUnit INTEGER (KIND = short), INTENT(IN) :: index !! time step !Local variables: TYPE (DateTime) :: time !!returned time INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: nDims !!number of dimensions INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: nAtts !!number of global attributes INTEGER (KIND = short) :: length !!length of time dimension INTEGER (KIND = short) :: idTime !!Id of the variable containing !!information on time ccordinate CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 100) :: variableName INTEGER (KIND = short) :: i !!loop index INTEGER :: slice (2) INTEGER :: timeSpan CHARACTER (LEN = 80) :: string CHARACTER (LEN = 19) :: str INTEGER, DIMENSION(NF90_MAX_VAR_DIMS) :: dimIDs !------------end of declaration------------------------------------------------ !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nDimensions = nDims, & nVariables = nVars, & nAttributes = nAtts ) CALL ncErrorHandler (ncStatus) !search for time variable DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN IF ( attribute(1:4) == 'time' ) THEN idTime = i EXIT END IF ELSE !standard_name is not defined: search for variable named 'time' !ncStatus = nf90_inq_varid (ncId, 'time', varid = i ) ncstatus = nf90_inquire_variable(ncId, varId = i, name = variableName) IF (LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'time' .OR. & LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'Time' .OR. & LEN_TRIM(variableName) == 5 .AND. & variableName(1:5) == 'Times') THEN !variable 'time' found idTime = i EXIT END IF END IF END DO !inquire time length ncStatus = nf90_inquire_variable(ncid, idTime, dimids = dimIDs) CALL ncErrorHandler (ncStatus) !ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(2), len = length) ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(1), len = length) CALL ncErrorHandler (ncStatus) !set next time step !IF (current < length) THEN ! next = current + 1 !ELSE ! next = length !END IF !compute date corresponding to next time step !slice(1) = 1 !slice(2) = next slice(1) = index slice(2) = 1 IF (DateTimeIsDefault(refTime)) THEN ncStatus = nf90_get_var (ncId, idTime, str , start = slice) CALL ncErrorHandler (ncStatus) string = str(1:10) // 'T' // str(12:19) // '+00:00' time = string ELSE ncStatus = nf90_get_var (ncId, idTime, timeSpan , start = slice) CALL ncErrorHandler (ncStatus) SELECT CASE (timeUnit) CASE ('minutes') timeSpan = timeSpan * minute CASE ('hours') timeSpan = timeSpan * hour CASE ('days') timeSpan = timeSpan * day CASE ('months') timeSpan = timeSpan * month END SELECT time = refTime + timeSpan END IF RETURN END FUNCTION TimeByIndex !============================================================================== !| Description: ! Parse units attribute of time variable ! Limits: ! string representing date and time must not contain blanks SUBROUTINE ParseTime & ! (ncId, ref_time, time_unit) USE StringManipulation, ONLY: & !Imported routines: StringCompact, StringToShort IMPLICIT NONE ! Arguments with intent(in): INTEGER (KIND = short), INTENT(IN) :: ncId !!NetCdf Id for the file ! Arguments with intent(out): TYPE(DateTime), INTENT(OUT) :: ref_time CHARACTER (LEN = 7), INTENT(OUT) :: time_unit ! local variables: INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: nDims !!number of dimensions INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: nAtts !!number of global attributes INTEGER (KIND = short) :: idTime !!Id of the variable containing !!information on time ccordinate CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 25) :: timeString CHARACTER (LEN = 80) :: unit CHARACTER (LEN = 10) :: since CHARACTER (LEN = 100) :: variableName LOGICAL :: error CHARACTER (LEN = 4) :: YYYY !! year CHARACTER (LEN = 2) :: MM !! month CHARACTER (LEN = 2) :: DD !! day CHARACTER (LEN = 2) :: hh !! hour CHARACTER (LEN = 2) :: min !! minute CHARACTER (LEN = 2) :: ss !! second CHARACTER (LEN = 6) :: tz !! time zone INTEGER (KIND = short) :: i !!loop index !------------end of declaration------------------------------------------------ !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nDimensions = nDims, & nVariables = nVars, & nAttributes = nAtts ) CALL ncErrorHandler (ncStatus) idTime = 0 DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN !standard_name is defined IF ( attribute(1:4) == 'time' ) THEN idTime = i EXIT END IF ELSE !standard_name is not defined: search for variable named 'time' or 'Times' !ncStatus = nf90_inq_varid (ncId, 'time', varid = i ) ncstatus = nf90_inquire_variable(ncId, varId = i, name = variableName) IF ( ( LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'time' ) .OR. & ( LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'Time' ) .OR. & ( LEN_TRIM(variableName) == 5 .AND. & variableName(1:5) == 'Times' ) ) THEN !variable found idTime = i EXIT END IF END IF END DO !get units attribute ncStatus = nf90_get_att (ncId, varid = idTime, name = 'units', & values = attribute) IF (ncStatus == nf90_noerr) THEN !units attribute found READ(attribute,*) unit i = INDEX (attribute, 'since', BACK = .TRUE.) IF (i > 0 ) THEN !attribute contains 'since' attribute = attribute ( i+6:LEN_TRIM(attribute) ) CALL ScanTimeStringAsString (attribute, YYYY, MM, DD, hh, min, ss, tz ) timeString = YYYY // '-' // MM // '-' // DD // 'T' // hh // ':' // min // ':' // ss // tz !set reference time. Convert to UTC ref_time = timeString ref_time = ToUTC (ref_time) ELSE !attribute is in the form like "day as %Y%m%d%h" ref_time = timeDefault END IF !set time_unit SELECT CASE (unit) CASE ('seconds', 'second', 'sec', 's') time_unit = 'seconds' CASE ('minutes', 'minute', 'min') time_unit = 'minutes' CASE ('hours', 'hour', 'hr', 'h') time_unit = 'hours' CASE ('days', 'day', 'd') time_unit = 'days' CASE ('months', 'month') time_unit = 'months' END SELECT ELSE ref_time = timeDefault END IF RETURN END SUBROUTINE !============================================================================== !| Description ! returns lower bounding time step of a given datetime SUBROUTINE SyncTime & ! (fileName, given, time) USE Units, ONLY: & ! Imported parameters: minute, hour, day, month USE StringManipulation, ONLY: & !Imported routines: ToString IMPLICIT NONE ! Arguments with intent(in): CHARACTER (LEN = *), INTENT (IN) :: fileName TYPE (DateTime), INTENT (IN) :: given ! Arguments with intent out TYPE (DateTime), INTENT (OUT) :: time !!returned time of the grid to sync to !Local variables: INTEGER (KIND = short) :: ncId !!NetCdf Id for the file TYPE (DateTime) :: refTime !!reference time to calculate time index CHARACTER (LEN = 30) :: timeUnit TYPE (DateTime) :: timeLower !!lower bound time TYPE (DateTime) :: timeUpper !!upper bound time !INTEGER (KIND = short) :: next !!index of next time step INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines INTEGER (KIND = short) :: nDims !!number of dimensions INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: nAtts !!number of global attributes INTEGER (KIND = short) :: length !!length of time dimension INTEGER (KIND = short) :: idTime !!Id of the variable containing !!information on time ccordinate CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 100) :: variableName INTEGER (KIND = short) :: i !!loop index INTEGER :: slice (2) INTEGER :: timeSpan CHARACTER (LEN = 80) :: string CHARACTER (LEN = 19) :: str INTEGER, DIMENSION(NF90_MAX_VAR_DIMS) :: dimIDs !------------end of declaration------------------------------------------------ !open file ncStatus = nf90_open (fileName, NF90_NOWRITE, ncId) IF (ncStatus /= nf90_noerr) THEN CALL Catch ('error', 'GridLib', & TRIM (nf90_strerror (ncStatus) )//': ', & code = ncIOError, argument = fileName ) ENDIF !Read time information CALL ParseTime (ncId, refTime, timeUnit) !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nDimensions = nDims, & nVariables = nVars, & nAttributes = nAtts ) CALL ncErrorHandler (ncStatus) !search for time variable DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN IF ( attribute(1:4) == 'time' ) THEN idTime = i EXIT END IF ELSE !standard_name is not defined: search for variable named 'time' !ncStatus = nf90_inq_varid (ncId, 'time', varid = i ) ncstatus = nf90_inquire_variable(ncId, varId = i, name = variableName) IF (LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'time' .OR. & LEN_TRIM(variableName) == 5 .AND. & variableName(1:5) == 'Times') THEN !variable 'time' found idTime = i EXIT END IF END IF END DO !inquire time length ncStatus = nf90_inquire_variable(ncid, idTime, dimids = dimIDs) CALL ncErrorHandler (ncStatus) !ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(2), len = length) ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(1), len = length) CALL ncErrorHandler (ncStatus) DO i = 1, length - 1 !set time step of lower bound slice(1) = i slice(2) = 1 !retrieve date and time ncStatus = nf90_get_var (ncId, idTime, timeSpan , start = slice) CALL ncErrorHandler (ncStatus) SELECT CASE (timeUnit) CASE ('minutes') timeSpan = timeSpan * minute CASE ('hours') timeSpan = timeSpan * hour CASE ('days') timeSpan = timeSpan * day CASE ('months') timeSpan = timeSpan * month END SELECT timeLower = refTime + timeSpan !set time step of upper bound slice(1) = i + 1 slice(2) = 1 !retrieve date and time ncStatus = nf90_get_var (ncId, idTime, timeSpan , start = slice) CALL ncErrorHandler (ncStatus) SELECT CASE (timeUnit) CASE ('minutes') timeSpan = timeSpan * minute CASE ('hours') timeSpan = timeSpan * hour CASE ('days') timeSpan = timeSpan * day CASE ('months') timeSpan = timeSpan * month END SELECT timeUpper = refTime + timeSpan IF ( given == timeLower ) THEN time = timeLower RETURN ELSE IF ( given == timeUpper ) THEN time = timeUpper RETURN ELSE IF ( given > timeLower .AND. given < timeUpper ) THEN time = timeLower RETURN END IF END DO !set next time step !IF (current < length) THEN ! next = current + 1 !ELSE ! next = length !END IF !compute date corresponding to next time step !slice(1) = 1 !slice(2) = next !slice(1) = next !slice(2) = 1 !IF (DateTimeIsDefault(refTime)) THEN ! ncStatus = nf90_get_var (ncId, idTime, str , start = slice) ! CALL ncErrorHandler (ncStatus) ! string = str(1:10) // 'T' // str(12:19) // '+00:00' ! time = string !ELSE ! ncStatus = nf90_get_var (ncId, idTime, timeSpan , start = slice) ! CALL ncErrorHandler (ncStatus) ! SELECT CASE (timeUnit) ! CASE ('minutes') ! timeSpan = timeSpan * minute ! CASE ('hours') ! timeSpan = timeSpan * hour ! CASE ('days') ! timeSpan = timeSpan * day ! CASE ('months') ! timeSpan = timeSpan * month ! END SELECT ! ! time = refTime + timeSpan ! !END IF RETURN END SUBROUTINE SyncTime !============================================================================== !| Description: ! get grid mapping for a floating point grid FUNCTION GetGridMappingFloat & ! ( layer ) & ! RESULT (grid_mapping) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_real), INTENT(INOUT) :: layer !Local variables TYPE (CRS) :: grid_mapping !------------end of declaration------------------------------------------------ grid_mapping = layer % grid_mapping END FUNCTION GetGridMappingFloat !============================================================================== !| Description: ! get grid mapping for a floating point grid FUNCTION GetGridMappingInteger & ! ( layer ) & ! RESULT (grid_mapping) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_integer), INTENT(INOUT) :: layer !Local variables TYPE (CRS) :: grid_mapping !------------end of declaration------------------------------------------------ grid_mapping = layer % grid_mapping END FUNCTION GetGridMappingInteger !============================================================================== !| Description: ! deallocate float grid SUBROUTINE GridDestroyFloat & ! (layer) IMPLICIT NONE ! Arguments with intent(inout): TYPE (grid_real), INTENT(INOUT) :: layer ! Local variables: INTEGER (KIND = short) :: ios !------------end of declaration------------------------------------------------ !deallocate data IF ( ALLOCATED (layer % mat) ) THEN DEALLOCATE ( layer % mat, STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory deallocating grid float ' , & code = memAllocError ) END IF END IF !deallocate coordinate reference system IF ( ALLOCATED ( layer % grid_mapping % param) ) THEN DEALLOCATE ( layer % grid_mapping % param, STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory deallocating grid float ' , & code = memAllocError ) END IF END IF IF ( ALLOCATED ( layer % grid_mapping % description) ) THEN DEALLOCATE ( layer % grid_mapping % description, STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory deallocating grid float ' , & code = memAllocError ) END IF END IF END SUBROUTINE GridDestroyFloat !============================================================================== !| Description: ! deallocate integer grid SUBROUTINE GridDestroyInteger & ! (layer) IMPLICIT NONE ! Arguments with intent(inout): TYPE (grid_integer), INTENT(INOUT) :: layer ! Local variables: INTEGER (KIND = short) :: ios !------------end of declaration------------------------------------------------ IF ( ALLOCATED (layer % mat) ) THEN DEALLOCATE ( layer % mat, STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory deallocating grid integer ' , & code = memAllocError ) END IF END IF !deallocate coordinate reference system IF ( ALLOCATED ( layer % grid_mapping % param) ) THEN DEALLOCATE ( layer % grid_mapping % param, STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory deallocating grid float ' , & code = memAllocError ) END IF END IF IF ( ALLOCATED ( layer % grid_mapping % description) ) THEN DEALLOCATE ( layer % grid_mapping % description, STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridLib', & 'memory deallocating grid float ' , & code = memAllocError ) END IF END IF END SUBROUTINE GridDestroyInteger !============================================================================== !| Description: ! transport matrix from netcdf format to grid_real SUBROUTINE SwapGridRealForward & ! (matIn, matOut, latlon) IMPLICIT NONE !arguments with intent in: REAL (KIND = float), INTENT (IN) :: matIn(:,:) !arguments with intent out: REAL (KIND = float), INTENT (OUT) :: matOut(:,:) !local variables: INTEGER :: i, j, idim, jdim, latlon !----------------------end of declaration-------------------------------------- idim = SIZE (matOut,1) jdim = SIZE (matOut,2) DO i = 1, idim DO j = 1, jdim IF (latlon == 1) THEN matOut (i,j) = matIn (j,idim - i + 1) ELSE matOut (i,j) = matIn (idim - i + 1,j) END IF END DO END DO END SUBROUTINE SwapGridRealForward !============================================================================== !| Description: ! transport matrix from netcdf format to grid_integer SUBROUTINE SwapGridIntegerForward & ! (matIn, matOut) IMPLICIT NONE !arguments with intent in: INTEGER (KIND = long), INTENT (IN) :: matIn(:,:) !arguments with intent out: INTEGER (KIND = long), INTENT (OUT) :: matOut(:,:) !local variables: INTEGER :: i,j,idim,jdim !----------------------end of declaration-------------------------------------- idim = SIZE (matOut,1) jdim = SIZE (matOut,2) DO i = 1, idim DO j = 1, jdim matOut (i,j) = matIn (j,idim - i + 1) END DO END DO END SUBROUTINE SwapGridIntegerForward !============================================================================== !| Description: ! transport matrix from grid_real to netcdf format SUBROUTINE SwapGridRealBack & ! (matIn, matOut) IMPLICIT NONE !arguments with intent in: REAL (KIND = float), INTENT (IN) :: matIn(:,:) !arguments with intent out: REAL (KIND = float), INTENT (OUT) :: matOut(:,:) !local variables: INTEGER :: i,j,idim,jdim !----------------------end of declaration-------------------------------------- idim = SIZE (matOut,1) jdim = SIZE (matOut,2) DO i = 1, idim DO j = 1, jdim matOut (i,j) = matIn (jdim - j + 1,i) END DO END DO END SUBROUTINE SwapGridRealBack !============================================================================== !| Description: ! transport matrix from grid_integer to netcdf format SUBROUTINE SwapGridIntegerBack & ! (matIn, matOut) IMPLICIT NONE !arguments with intent in: INTEGER (KIND = long), INTENT (IN) :: matIn(:,:) !arguments with intent out: INTEGER (KIND = long), INTENT (OUT) :: matOut(:,:) !local variables: INTEGER :: i,j,idim,jdim !----------------------end of declaration-------------------------------------- idim = SIZE (matOut,1) jdim = SIZE (matOut,2) DO i = 1, idim DO j = 1, jdim matOut (i,j) = matIn (jdim - j + 1,i) END DO END DO END SUBROUTINE SwapGridIntegerBack !============================================================================== !| Description: ! Handle errors from netcdf related operation SUBROUTINE ncErrorHandler & ! (errcode) IMPLICIT NONE ! Local scalars: INTEGER (KIND = short), INTENT (IN) :: errcode !------------end of declaration------------------------------------------------ IF (errcode /= nf90_noerr) THEN CALL Catch ('error', 'GridLib', & TRIM (nf90_strerror (errcode) ), & code = ncIOError ) ENDIF END SUBROUTINE ncErrorHandler !! Description: !! given filename of a multidimensional net-cdf file !! the GetDtGrid function returns the temporal resolution (seconds) !! assuming that it is regular !! If option checkRegular is true the function check that !! temporal resolution is regular FUNCTION GetDtGrid & ! (filename, checkRegular) & ! RESULT (dt) USE Units, ONLY: & ! Imported parameters: minute, hour, day, month USE StringManipulation, ONLY: & ! imported routines: ToString IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: filename LOGICAL, OPTIONAL, INTENT(IN) :: checkRegular !Local declarations: INTEGER (KIND = long) :: dt INTEGER (KIND = short) :: ncStatus !!error code returned by NetCDF routines INTEGER (KIND = short) :: ncId !!NetCdf Id for the file INTEGER (KIND = short) :: nDims !!number of dimensions INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: nAtts !!number of global attributes INTEGER (KIND = short) :: timeDimId !!id of time dimension INTEGER (KIND = short) :: length !!length of time dimension INTEGER (KIND = short) :: idTime !!Id of the variable containing !!information on time coordinate CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 25) :: timeString CHARACTER (LEN = 100) :: variableName TYPE(DateTime) :: ref_time TYPE(DateTime) :: date1, date2 CHARACTER (LEN = 7) :: time_unit INTEGER :: slice (1) INTEGER :: slice2 (2) INTEGER :: time1, time2 INTEGER :: timeSpan INTEGER (KIND = short) :: i INTEGER, DIMENSION(NF90_MAX_VAR_DIMS) :: dimIDs CHARACTER (LEN = 19) :: str, str1, str2 !------------end of declaration------------------------------------------------ !open net-cdf file with read-only access ncStatus = nf90_open (fileName, NF90_NOWRITE, ncId) IF (ncStatus /= nf90_noerr) THEN CALL Catch ('error', 'GridLib', & TRIM (nf90_strerror (ncStatus) )//': ', & code = ncIOError, argument = fileName ) ENDIF !retrieve time unit CALL ParseTime (ncId, ref_time, time_unit) !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nDimensions = nDims, & nVariables = nVars, & nAttributes = nAtts ) CALL ncErrorHandler (ncStatus) !search for time variable DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN !standard_name is defined IF ( attribute(1:4) == 'time' ) THEN idTime = i EXIT END IF ELSE !standard_name is not defined: search for variable named 'time' !ncStatus = nf90_inq_varid (ncId, 'time', varid = i ) ncstatus = nf90_inquire_variable(ncId, varId = i, name = variableName) IF (LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'time' .OR. & LEN_TRIM(variableName) == 5 .AND. & variableName(1:5) == 'Times' ) THEN !variable 'time' found idTime = i EXIT END IF END IF END DO !retrieve time length length = GetTimeSteps (ncId) !compute dt IF (DateTimeIsDefault(ref_time)) THEN slice2(1) = 1 slice2(2) = 1 ncStatus = nf90_get_var (ncId, idTime, str , start = slice2) CALL ncErrorHandler (ncStatus) !build datetime string from format used in netcdf file i.e 2007-10-11_00:00:00 timeString = str(1:10) // 'T' // str(12:19) // '+00:00' date1 = timeString slice2(1) = 1 slice2(2) = 2 ncStatus = nf90_get_var (ncId, idTime, str , start = slice2) CALL ncErrorHandler (ncStatus) !build datetime string from format used in netcdf file i.e 2007-10-11_00:00:00 timeString = str(1:10) // 'T' // str(12:19) // '+00:00' date2 = timeString dt = date2 - date1 ELSE slice(1) = 1 ncStatus = nf90_get_var (ncId, idTime, time1 , start = slice) CALL ncErrorHandler (ncStatus) slice(1) = 2 ncStatus = nf90_get_var (ncId, idTime, time2 , start = slice) CALL ncErrorHandler (ncStatus) dt = time2 - time1 !convert in seconds SELECT CASE (time_unit) CASE ('minutes') dt = dt * minute CASE ('hours') dt = dt * hour CASE ('days') dt = dt * day CASE ('months') dt = dt * month END SELECT END IF !Check if dt is regular IF (PRESENT (checkRegular) ) THEN IF (checkRegular) THEN DO i = 1, length - 1 IF (DateTimeIsDefault(ref_time)) THEN slice2(1) = 1 slice2(2) = i ncStatus = nf90_get_var (ncId, idTime, str1 , start = slice2) CALL ncErrorHandler (ncStatus) slice2(2) = i + 1 ncStatus = nf90_get_var (ncId, idTime, str2 , start = slice2) CALL ncErrorHandler (ncStatus) timeString = str1(1:10) // 'T' // str1(12:19) // '+00:00' date1 = timeString timeString = str2(1:10) // 'T' // str2(12:19) // '+00:00' date2 = timeString timeSpan = date2 - date1 ELSE slice(1) = i ncStatus = nf90_get_var (ncId, idTime, time1 , start = slice) CALL ncErrorHandler (ncStatus) slice(1) = i + 1 ncStatus = nf90_get_var (ncId, idTime, time2 , start = slice) CALL ncErrorHandler (ncStatus) IF (DateTimeIsDefault(ref_time)) THEN timeString = ToString (time1) timeString = timeString (1:4) // '-' // & timeString (5:6) // '-' // & timeString (7:8) // 'T' // & timeString (9:10) // ':00:00+00:00' date1 = timeString timeString = ToString (time2) timeString = timeString (1:4) // '-' // & timeString (5:6) // '-' // & timeString (7:8) // 'T' // & timeString (9:10) // ':00:00+00:00' date2 = timeString timeSpan = date2 - date1 ELSE timeSpan = time2 - time1 !convert in seconds SELECT CASE (time_unit) CASE ('minutes') timeSpan = timeSpan * minute CASE ('hours') timeSpan = timeSpan * hour CASE ('days') timeSpan = timeSpan * day CASE ('months') timeSpan = timeSpan * month END SELECT END IF END IF IF (timeSpan /= dt ) THEN CALL Catch ('error', 'GridLib', & 'time not regular in multidimensional grid') END IF END DO END IF END IF END FUNCTION GetDtGrid !============================================================================== !| Description: ! given filename of a multidimensional net-cdf file ! the GetFirstDate function returns the date and time ! of first grid FUNCTION GetFirstDate & ! (filename) & ! RESULT (time) USE Chronos, ONLY: & !Imported type definitions: DateTime USE Units, ONLY: & ! Imported parameters: minute, hour, day, month IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: filename !Local declarations: TYPE (DateTime) :: time INTEGER (KIND = short) :: ncStatus !!error code returned by NetCDF routines INTEGER (KIND = short) :: ncId !!NetCdf Id for the file TYPE(DateTime) :: ref_time CHARACTER (LEN = 7) :: time_unit INTEGER (KIND = short) :: idTime !!Id of the variable containing !!information on time coordinate INTEGER (KIND = short) :: nDims !!number of dimensions INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: nAtts !!number of global attributes CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 100) :: variableName INTEGER :: slice (1) INTEGER :: slice2 (2) INTEGER :: time1 INTEGER :: dt INTEGER (KIND = short) :: i CHARACTER (LEN = 19) :: str !------------end of declaration------------------------------------------------ !open net-cdf file with read-only access ncStatus = nf90_open (fileName, NF90_NOWRITE, ncId) IF (ncStatus /= nf90_noerr) THEN CALL Catch ('error', 'GridLib', & TRIM (nf90_strerror (ncStatus) )//': ', & code = ncIOError, argument = fileName ) ENDIF !retrieve time unit CALL ParseTime (ncId, ref_time, time_unit) !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nDimensions = nDims, & nVariables = nVars, & nAttributes = nAtts ) CALL ncErrorHandler (ncStatus) !search for time variable DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN !standard_name is defined IF ( attribute(1:4) == 'time' ) THEN idTime = i EXIT END IF ELSE !standard_name is not defined: search for variable named 'time' !ncStatus = nf90_inq_varid (ncId, 'time', varid = i ) ncstatus = nf90_inquire_variable(ncId, varId = i, name = variableName) IF (LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'time' .OR. & LEN_TRIM(variableName) == 5 .AND. & variableName(1:5) == 'Times' ) THEN !variable 'time' found idTime = i EXIT END IF END IF END DO IF (DateTimeIsDefault(ref_time)) THEN slice2(1) = 1 slice2(2) = 1 ncStatus = nf90_get_var (ncId, idTime, str , start = slice2) CALL ncErrorHandler (ncStatus) !build datetime string from format used in netcdf file i.e 2007-10-11_00:00:00 timeString = str(1:10) // 'T' // str(12:19) // '+00:00' time = timeString ELSE !retrieve timespan of first grid slice(1) = 1 ncStatus = nf90_get_var (ncId, idTime, time1 , start = slice) CALL ncErrorHandler (ncStatus) dt = time1 !convert in seconds SELECT CASE (time_unit) CASE ('minutes') dt = dt * minute CASE ('hours') dt = dt * hour CASE ('days') dt = dt * day CASE ('months') dt = dt * month END SELECT !compute time time = ref_time + dt END IF RETURN END FUNCTION GetFirstDate !============================================================================== !| Description: ! given filename of a multidimensional net-cdf file ! the GetTimeStepsFromFile function returns the number of time steps FUNCTION GetTimeStepsFromFile & ! (filename) & ! RESULT (steps) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: filename !Local declarations: INTEGER (KIND = short) :: steps INTEGER (KIND = short) :: ncStatus !!error code returned by NetCDF routines INTEGER (KIND = short) :: ncId !!NetCdf Id for the file INTEGER (KIND = short) :: nDims !!number of dimensions INTEGER (KIND = short) :: nDimsVar !!number of dimensions of a variable INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: nAtts !!number of global attributes INTEGER (KIND = short) :: idTime !!Id of the variable containing !!information on time coordinate INTEGER (KIND = short) :: dimId !dimension id CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 100) :: variableName INTEGER (KIND = short) :: i INTEGER, DIMENSION(NF90_MAX_VAR_DIMS) :: dimIDs !------------end of declaration------------------------------------------------ !open net-cdf file with read-only access ncStatus = nf90_open (fileName, NF90_NOWRITE, ncId) IF (ncStatus /= nf90_noerr) THEN CALL Catch ('error', 'GridLib', & TRIM (nf90_strerror (ncStatus) )//': ', & code = ncIOError, argument = fileName ) ENDIF !search for "Time" dimension ncStatus = nf90_inq_dimid(ncid, "Time", dimId) IF (ncStatus == nf90_noerr) THEN !Time dimension found ncStatus = nf90_inquire_dimension (ncId, dimId, len = steps) RETURN !no need to go on END IF !search for "time" dimension ncStatus = nf90_inq_dimid(ncid, "time", dimId) IF (ncStatus /= nf90_noerr) THEN !time dimension found ncStatus = nf90_inquire_dimension (ncId, dimId, len = steps) RETURN !no need to go on END IF !If dimension has a different name, analyse variable !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nDimensions = nDims, & nVariables = nVars, & nAttributes = nAtts ) CALL ncErrorHandler (ncStatus) !search for time variable DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN !standard_name is defined IF ( attribute(1:4) == 'time' ) THEN idTime = i EXIT END IF ELSE !standard_name is not defined: search for variable named 'time' !ncStatus = nf90_inq_varid (ncId, 'time', varid = i ) ncstatus = nf90_inquire_variable(ncId, varId = i, name = variableName) IF (LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'time' .OR. & LEN_TRIM(variableName) == 5 .AND. & variableName(1:5) == 'Times' ) THEN !variable 'time' found idTime = i EXIT END IF END IF END DO !retrieve dimensions of variable ncStatus = nf90_inquire_variable (ncId, idTime, ndims = nDimsVar) CALL ncErrorHandler (ncStatus) !retrieve dimension id of variable ncStatus = nf90_inquire_variable(ncid, idTime, dimids = dimIDs) CALL ncErrorHandler (ncStatus) !inquire time length IF (nDimsVar == 1) THEN ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(1), len = steps) CALL ncErrorHandler (ncStatus) ELSE IF (nDimsVar == 2) THEN ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(2), len = steps) CALL ncErrorHandler (ncStatus) END IF RETURN END FUNCTION GetTimeStepsFromFile !============================================================================== !| Description: ! given ncId of a multidimensional net-cdf file ! the GetTimeStepsFromNCid function returns the number of time steps FUNCTION GetTimeStepsFromNCid & ! (ncId) & ! RESULT (steps) IMPLICIT NONE !Arguments with intent(in): INTEGER (KIND = short), INTENT(IN) :: ncId !!NetCdf Id for the file !Local declarations: INTEGER (KIND = short) :: steps INTEGER (KIND = short) :: ncStatus !!error code returned by NetCDF routines INTEGER (KIND = short) :: nDims !!number of dimensions INTEGER (KIND = short) :: nDimsVar !!number of dimensions of a variable INTEGER (KIND = short) :: nVars !!number of variables INTEGER (KIND = short) :: nAtts !!number of global attributes INTEGER (KIND = short) :: idTime !!Id of the variable containing !!information on time coordinate INTEGER (KIND = short) :: dimId !dimension id CHARACTER (LEN = 80) :: attribute CHARACTER (LEN = 100) :: variableName INTEGER (KIND = short) :: i INTEGER, DIMENSION(NF90_MAX_VAR_DIMS) :: dimIDs !------------end of declaration------------------------------------------------ !search for "Time" dimension ncStatus = nf90_inq_dimid(ncid, "Time", dimId) IF (ncStatus == nf90_noerr) THEN !Time dimension found ncStatus = nf90_inquire_dimension (ncId, dimId, len = steps) RETURN !no need to go on END IF !search for "time" dimension ncStatus = nf90_inq_dimid(ncid, "time", dimId) IF (ncStatus /= nf90_noerr) THEN !time dimension found ncStatus = nf90_inquire_dimension (ncId, dimId, len = steps) RETURN !no need to go on END IF !If dimension has a different name, analyse variable !inquire dataset to retrieve number of dimensions, variables !and global attributes ncStatus = nf90_inquire(ncId, nDimensions = nDims, & nVariables = nVars, & nAttributes = nAtts ) CALL ncErrorHandler (ncStatus) !search for time variable DO i = 1, nVars attribute = '' ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', & values = attribute) IF (ncStatus == nf90_noerr) THEN !standard_name is defined IF ( attribute(1:4) == 'time' ) THEN idTime = i EXIT END IF ELSE !standard_name is not defined: search for variable named 'time' !ncStatus = nf90_inq_varid (ncId, 'time', varid = i ) ncstatus = nf90_inquire_variable(ncId, varId = i, name = variableName) IF (LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'time' .OR. & LEN_TRIM(variableName) == 4 .AND. & variableName(1:4) == 'Time' .OR. & LEN_TRIM(variableName) == 5 .AND. & variableName(1:5) == 'Times' ) THEN !variable 'time' found idTime = i EXIT END IF END IF END DO !retrieve dimensions of variable ncStatus = nf90_inquire_variable (ncId, idTime, ndims = nDimsVar) CALL ncErrorHandler (ncStatus) !retrieve dimension id of variable ncStatus = nf90_inquire_variable(ncid, idTime, dimids = dimIDs) CALL ncErrorHandler (ncStatus) !inquire time length IF (nDimsVar == 1) THEN ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(1), len = steps) CALL ncErrorHandler (ncStatus) ELSE IF (nDimsVar == 2) THEN ncStatus = nf90_inquire_dimension (ncId, dimid = dimIDs(2), len = steps) CALL ncErrorHandler (ncStatus) END IF RETURN END FUNCTION GetTimeStepsFromNCid !============================================================================== !| Description: ! scan a string to extract time information as character strings: ! year, month, day, hour, minute, second, timezone ! supported formats: ! ! * `1900-01-01T01:00:00+00:00` ! * `1900-1-1 01:00:00` ! * `1900-1-1 1:0:0` ! SUBROUTINE ScanTimeStringAsString & ! (string, YYYY, MM, DD, hh, min, ss, tz ) USE StringManipulation, ONLY: & !Imported routines: StringTokenize, StringCompact, & StringToShort, ToString IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT (IN) :: string CHARACTER (LEN = *), INTENT (OUT) :: YYYY !! year CHARACTER (LEN = *), INTENT (OUT) :: MM !! month CHARACTER (LEN = *), INTENT (OUT) :: DD !! day CHARACTER (LEN = *), INTENT (OUT) :: hh !! hour CHARACTER (LEN = *), INTENT (OUT) :: min !! minute CHARACTER (LEN = *), INTENT (OUT) :: ss !! second CHARACTER (LEN = *), INTENT (OUT) :: tz !!time zone !Local declarations: CHARACTER (LEN = 2) :: tzhh !! time zone hour CHARACTER (LEN = 2) :: tzmin !! time zone minute INTEGER (KIND = short) :: i, iplus, iminus INTEGER (KIND = short) :: year, month, day, hour, minute, second, tzhour, tzminute CHARACTER (len=100), POINTER :: args (:) CHARACTER (LEN = 20) :: time, tzstring INTEGER (KIND = short) :: nargs character(len=1024) :: format_string !------------end of declaration------------------------------------------------ !search for year, month, and day CALL StringTokenize (string = string, delims = '-', args = args, nargs = nargs) year = StringToShort ( StringCompact (args (1)) ) month = StringToShort ( StringCompact (args (2)) ) !search for 'T' after day i = INDEX (args (3), 'T', BACK = .TRUE.) IF (i > 0) THEN day = StringToShort ( args (3) (1 : i-1) ) ELSE day = StringToShort ( StringCompact (args (3)) ) END IF !search for time: hour, minute, second !check wether T is present IF ( INDEX (string, 'T') > 0) THEN i = INDEX (string, 'T') ELSE IF ( INDEX (string, ' ') > 0) THEN i = INDEX (string, ' ') END IF time = string ( i+1 : len_trim (string) ) CALL StringTokenize (string = time, delims = ':', args = args, nargs = nargs) hour = StringToShort ( StringCompact (args (1)) ) minute = StringToShort ( StringCompact (args (2)) ) iplus = INDEX (args(3), '+') iminus = INDEX (args(3), '-') format_string = "(I2.2)" IF ( iplus > 0) THEN second = StringToShort ( StringCompact (args (3) (1 : iplus - 1) ) ) CALL StringTokenize (string = time, delims = '+', args = args, nargs = nargs) tzstring = args (2) CALL StringTokenize (string = tzstring, delims = ':', args = args, nargs = nargs) tzhour = StringToShort ( StringCompact (args (1)) ) tzminute = StringToShort ( StringCompact (args (2)) ) tzhh = ToString (tzhour, format_string) tzmin = ToString (tzminute, format_string) tz = '+' // tzhh // ':' // tzmin ELSE IF (iminus > 0 ) THEN second = StringToShort ( StringCompact (args (3) (1 : iminus - 1) ) ) CALL StringTokenize (string = time, delims = '-', args = args, nargs = nargs) tzstring = args (2) CALL StringTokenize (string = tzstring, delims = ':', args = args, nargs = nargs) tzhour = StringToShort ( StringCompact (args (1)) ) tzminute = StringToShort ( StringCompact (args (2)) ) tzhh = ToString (tzhour, format_string) tzmin = ToString (tzminute, format_string) tz = '-' // tzhh // ':' // tzmin ELSE second = StringToShort ( StringCompact (args (3) ) ) !time zone is not given, set to default tzhh = '00' tzmin = '00' tz = '+' // tzhh // ':' // tzmin END IF !convert short integer to string YYYY = ToString (year) MM = ToString (month, format_string) DD = ToString (day, format_string) hh = ToString (hour, format_string) min = ToString (minute, format_string) ss = ToString (second, format_string) RETURN END SUBROUTINE ScanTimeStringAsString END MODULE GridLib